diff --git a/holodeck/anisotropy.py b/holodeck/anisotropy.py index 6dc340c1..b6725f99 100644 --- a/holodeck/anisotropy.py +++ b/holodeck/anisotropy.py @@ -53,7 +53,7 @@ def healpix_map(hc_ss, hc_bg, nside=NSIDE, seed=None, ret_seed=False): np.random.seed(seed) # spread background evenly across pixels in moll_hc - moll_hc = np.ones((nfreqs,nreals,npix)) * hc_bg[:,:,np.newaxis]**2/(npix/area) # (frequency, realization, pixel) + moll_hc = np.ones((nfreqs,nreals,npix)) * hc_bg[:,:,np.newaxis]**2/(npix*area) # (frequency, realization, pixel) # choose random pixels to place the single sources pix_ss = np.random.randint(0, npix-1, size=nfreqs*nreals*nloudest).reshape(nfreqs, nreals, nloudest) diff --git a/holodeck/cyutils.pyx b/holodeck/cyutils.pyx index 08f3dc08..9dc3f448 100644 --- a/holodeck/cyutils.pyx +++ b/holodeck/cyutils.pyx @@ -23,7 +23,7 @@ cimport scipy.special.cython_special as sp_special from libc.stdio cimport printf from libc.stdlib cimport malloc, free, qsort # make sure to use c-native math functions instead of python/numpy -from libc.math cimport pow, sqrt, abs, M_PI, NAN +from libc.math cimport pow, sqrt, abs, M_PI, NAN, cos, sin from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer from numpy.random cimport bitgen_t @@ -1815,4 +1815,285 @@ cdef double[:, :] _interp_2d( else: ynew[ii, nn] = NAN - return ynew \ No newline at end of file + return ynew + + + + + + + + +# ================================================================================================== +# ==== DetStats Functions ==== +# ================================================================================================== + + +def gamma_of_rho_interp(rho, rsort, rho_interp_grid, gamma_interp_grid): + """ + rho : 1Darray of scalars + SNR of single sources, in flat array + rsort : 1Darray + order of flat rho values smallest to largest + rho_interp_grid : 1Darray + rho values corresponding to each gamma + gamma_interp_grid : 1Darray + gamma values corresponding to each rho + + """ + # pass in the interp grid + cdef np.ndarray[np.double_t, ndim=1] gamma = np.zeros(rho.shape) + + _gamma_of_rho_interp(rho, rsort, rho_interp_grid, gamma_interp_grid, gamma) + + return gamma + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.nonecheck(False) +@cython.cdivision(True) +cdef int _gamma_of_rho_interp( + double[:] rho, long[:] rsort, + double[:] rho_interp_grid, double[:] gamma_interp_grid, + # output + double[:] gamma + ): + """ Find gamma of rho by interpolation over rho and gamma grids. + """ + + cdef int n_rho = rho.size + cdef int n_interp = rho_interp_grid.size + cdef int ii, kk, rr + ii = 0 # get rho in order using rho[rsort[ii]] + + for kk in range(n_rho): + rr = rsort[kk] # index of next largest rho, equiv to rev in redz calculation + # print('kk =',kk,' rr =', rr, 'rho[rr] =', rho[rr]) + # get to the right index of the interpolation-grid + while (rho_interp_grid[ii+1] < rho[rr]) and (ii < n_interp -2): + ii += 1 + # print('ii =',ii, ' rho_interp[ii] =', rho_interp_grid[ii], ' rho_interp[ii+1] =', rho_interp_grid[ii+1]) + # interpolate + gamma[rr] = interp_at_index(ii, rho[rr], rho_interp_grid, gamma_interp_grid) + # print('rho =', rho[rr], ' gamma =', gamma[rr], '\n')\ + + return 0 + + +def snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): + """ Calculate single source SNR + + + Parameters + ---------- + amp : (F,R,L) NDarray + Dimensionless strain amplitude of loudest single sources + F_iplus : (P,F,S,L) NDarray + Antenna pattern function for each pulsar. + F_icross : (P,F,S,L) NDarray + Antenna pattern function for each pulsar. + iotas : (F,S,L) NDarray + Inclination, used to calculate: + a_pol = 1 + np.cos(iotas) **2 + b_pol = -2 * np.cos(iotas) + dur : scalar + Duration of observations. + Phi_0 : (F,S,L) NDarray + Initial GW phase + S_i : (P,F,R,L) NDarray + Total noise of each pulsar wrt detection of each single source, in s^3 + freqs : (F,) 1Darray + Observed frequency bin centers. + + Returns + ------- + snr_ss : (F,R,S,L) NDarray + SNR from the whole PTA for each single source with + each realized sky position (S) and realized strain (R) + + """ + nfreqs, nreals, nloudest = amp.shape[0], amp.shape[1], amp.shape[2] + npsrs, nskies = F_iplus.shape[0], F_iplus.shape[2] + cdef np.ndarray[np.double_t, ndim=4] snr_ss = np.zeros((nfreqs, nreals, nskies, nloudest)) + _snr_ss( + amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs, + npsrs, nfreqs, nreals, nskies, nloudest, + snr_ss) + return snr_ss + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.nonecheck(False) +@cython.cdivision(True) +cdef int _snr_ss( + double[:,:,:] amp, + double[:,:,:,:] F_iplus, + double[:,:,:,:] F_icross, + double[:,:,:] iotas, + double dur, + double[:,:,:] Phi_0, + double[:,:,:,:] S_i, + double[:] freqs, + long npsrs, long nfreqs, long nreals, long nskies, long nloudest, + # output + double[:,:,:,:] snr_ss + ): + """ + + Parameters + ---------- + amp : (F,R,L) NDarray + Dimensionless strain amplitude of loudest single sources + F_iplus : (P,F,S,L) NDarray + Antenna pattern function for each pulsar. + F_icross : (P,F,S,L) NDarray + Antenna pattern function for each pulsar. + iotas : (F,S,L) NDarray + Inclination, used to calculate: + a_pol = 1 + np.cos(iotas) **2 + b_pol = -2 * np.cos(iotas) + dur : scalar + Duration of observations. + Phi_0 : (F,S,L) NDarray + Initial GW phase + S_i : (P,F,R,L) NDarray + Total noise of each pulsar wrt detection of each single source, in s^3 + freqs : (F,) 1Darray + Observed frequency bin centers. + snr_ss : (F,R,S,L) NDarray + Pointer to single source SNR array, to be calculated. + + NOTE: This may be improved by moving some of the math outside the function. + I.e., passing in sin/cos of NDarrays to be used. + """ + + cdef int pp, ff, rr, ss, ll + cdef float a_pol, b_pol, Phi_T, pta_snr_sq, coef, term1, term2, term3 + # print('npsrs %d, nfreqs %d, nreals %d, nskies %d, nloudest %d' % (npsrs, nfreqs, nreals, nskies, nloudest)) + + for ff in range(nfreqs): + for ss in range(nskies): + for ll in range(nloudest): + a_pol = 1 + pow(cos(iotas[ff,ss,ll]), 2.0) + b_pol = -2 * cos(iotas[ff,ss,ll]) + Phi_T = 2 * M_PI * freqs[ff] * dur + Phi_0[ff,ss,ll] + for rr in range(nreals): + pta_snr_sq = 0 + for pp in range(npsrs): + # calculate coefficient depending on + # function of amp, S_i, and freqs + coef = pow(amp[ff,rr,ll], 2.0) / (S_i[pp,ff,rr,ll] * 8 * pow(M_PI * freqs[ff], 3.0)) + + # calculate terms that depend on p, f, s, and l + # functions of F_iplus, F_icross, a_pol, b_pol, Phi_0, and Phi_T + term1 = ( + pow(a_pol * F_iplus[pp,ff,ss,ll], 2.0) + * (Phi_T * (1.0 + 2.0 * pow(sin(Phi_0[ff,ss,ll]), 2.0)) + + cos(Phi_T) * (-1.0 * sin(Phi_T) + 4.0 * cos(Phi_0[ff,ss,ll])) + - 4.0 * sin(Phi_0[ff,ss,ll]) + ) + ) + term2 = ( + pow(b_pol * F_icross[pp,ff,ss,ll], 2.0) + * (Phi_T * (1.0 + 2.0 * pow(cos(Phi_0[ff,ss,ll]), 2.0)) + + sin(Phi_T) * cos(Phi_T) - 4.0 * cos(Phi_0[ff,ss,ll]) + ) + ) + term3 = ( + -2.0 * a_pol * b_pol * F_iplus[pp,ff,ss,ll] * F_icross[pp,ff,ss,ll] + * (2.0 * Phi_T * sin(Phi_T) *cos(Phi_0[ff,ss,ll]) + + sin(Phi_T) * (sin(Phi_T) - 2.0 * sin(Phi_0[ff,ss,ll]) + + 2.0 * cos(Phi_T) * cos(Phi_0[ff,ss,ll]) + - 2.0 * cos(Phi_0[ff,ss,ll]) + ) + ) + ) + pta_snr_sq += coef*(term1 + term2 + term3) # sum snr^2 of all pulsars for a single source + + # set snr for a single source, using sum from all pulsars + snr_ss[ff,rr,ss,ll] = sqrt(pta_snr_sq) + + +def Sh_rest(hc_ss, hc_bg, freqs, nexcl): + """ + Calculate the noise from all the single sources except the source in question + and the next N_excl loudest sources. + + Parameters + ---------- + hc_ss : (F,R,L) NDarray + Characteristic strain from all loud single sources. + hc_bg : (F,R) NDarray + Characteristic strain from all but loudest source at each frequency. + freqs : (F,) 1Darray + Frequency bin centers. + nexcl : int + Number of loudest single sources to exclude from hc_rest noise, in addition + to the source in question. + + Returns + ------- + Sh_rest : (F,R,L) NDarray of scalars + The noise in a single pulsar from other GW sources for detecting each single source. + + """ + + nfreqs, nreals, nloudest = hc_ss.shape[0], hc_ss.shape[1], hc_ss.shape[2] + cdef np.ndarray[np.double_t, ndim=3] Sh_rest = np.zeros((nfreqs, nreals, nloudest)) + _Sh_rest(hc_ss, hc_bg, freqs, nexcl, nreals, nfreqs, nloudest, Sh_rest) + return Sh_rest + + +cdef void _Sh_rest( + double[:,:,:] hc_ss, double[:,:,] hc_bg, double[:] freqs, long nexcl, + long nreals, long nfreqs, long nloudest, + double[:,:,:] Sh_rest): + """ + Calculate the noise from all the single sources except the source in question + and the next N_excl loudest sources. + + Parameters + ---------- + hc_ss : (F,R,L) NDarray + Characteristic strain from all loud single sources. + hc_bg : (F,R) NDarray + Characteristic strain from all but loudest source at each frequency. + freqs : (F,) 1Darray + Frequency bin centers. + nexcl : int + Number of loudest single sources to exclude from hc_rest noise, in addition + to the source in question. + + Returns + ------- + void + Sh_rest : (F,R,L) NDarray of scalars, updated via memory address + + + Sh = hc^2 / (f^3 12 pi^2) + + """ + cdef int ff, rr, ll, count + cdef double Sh_ss, Sh_bg + + for ff in range(nfreqs): + freq = freqs[ff] + for rr in range(nreals): + for ll in range(nloudest): # calculating for the llth loduest + Sh_ss = 0 + count = 0 + for ii in range(nloudest): # adding other loudest with index ii + if (ii != ll): # check it's not our current source + # if current is in top N_excl, must be (N_excl+1)th or above + # if current is >= top N_excl, must be (N_excl)th or above + if ((ll < nexcl) and (ii > nexcl)) or ((ll >= nexcl) and (ii > nexcl-1)): + Sh_ss += pow(hc_ss[ff,rr,ii], 2.0) / pow(freq, 3.0) / (12*pow(M_PI, 2.0)) + count += 1 + Sh_bg = pow(hc_bg[ff,rr], 2.0) / pow(freq, 3.0) / (12*pow(M_PI, 2.0)) + Sh_rest[ff,rr,ll] = Sh_rest[ff,rr,ll] + Sh_ss + Sh_bg + if count != (nloudest - nexcl - 1): + err = (f"ERROR in calculate Sh_rest! count of sources={count}") + print(err) + + + diff --git a/holodeck/detstats.py b/holodeck/detstats.py index 78e40543..50191650 100644 --- a/holodeck/detstats.py +++ b/holodeck/detstats.py @@ -12,12 +12,13 @@ import matplotlib.pyplot as plt import os from datetime import datetime +import warnings import holodeck as holo -from holodeck import utils, cosmo, log, plot, sam_cython +from holodeck import utils, cosmo, log, plot, anisotropy from holodeck.constants import MPC, YR -from holodeck.sams import cyutils as sam_cyutils +from holodeck import cyutils import hasasia.sensitivity as hsen import hasasia.sim as hsim @@ -183,7 +184,7 @@ def _white_noise(delta_t, sigma_i): ######################## Power Spectral Density ######################## -def _power_spectral_density(hc_bg, freqs): +def _power_spectral_density(hc_bg, freqs, reshape_freqs=True): """ Calculate the spectral density S_h(f_k) ~ S_h0(f_k) at the kth frequency Parameters @@ -202,8 +203,10 @@ def _power_spectral_density(hc_bg, freqs): Follows Eq. (25) of Rosado et al. 2015 """ + if reshape_freqs: + freqs = freqs[:,np.newaxis] - S_h = hc_bg**2 / (12 * np.pi**2 * freqs[:,np.newaxis]**3) + S_h = hc_bg**2 / (12 * np.pi**2 * freqs**3) return S_h @@ -488,8 +491,9 @@ def detect_bg(thetas, phis, sigmas, fobs, cad, hc_bg, alpha_0=0.001, ret = False -def detect_bg_pta(pulsars, fobs, hc_bg, hc_ss, alpha_0=0.001, ret_snr = False, - red_amp=None, red_gamma=None, ss_noise=True): +def detect_bg_pta(pulsars, fobs, hc_bg, hc_ss=None, custom_noise=None, + alpha_0=0.001, ret_snr = False, + red_amp=None, red_gamma=None, ss_noise=False): """ Calculate the background detection probability, and all the intermediary steps from a list of hasasia.Pulsar objects. @@ -538,18 +542,25 @@ def detect_bg_pta(pulsars, fobs, hc_bg, hc_ss, alpha_0=0.001, ret_snr = False, Sh_bg = _power_spectral_density(hc_bg[:], fobs) Sh0_bg = Sh_bg # note this refers to same object, not a copy - # calculate white noise - noise = _white_noise(cad, sigmas)[:,np.newaxis] # P,1 + # noise spectral density + if custom_noise is not None: + if custom_noise.shape != (len(pulsars), len(fobs), len(hc_bg[0])): + err = f"{custom_noise.shape=}, must be shape (P,F,R)=({len(pulsars)}, {len(fobs)}, {len(hc_bg[0])})" + raise ValueError(err) + noise = custom_noise + else: + # calculate white noise + noise = _white_noise(cad, sigmas)[:,np.newaxis] # P,1 - # add red noise - if (red_amp is not None) and (red_gamma is not None): - red_noise = _red_noise(red_amp, red_gamma, fobs)[np.newaxis,:] # (1,F,) - noise = noise + red_noise # (P,F,) + # add red noise + if (red_amp is not None) and (red_gamma is not None): + red_noise = _red_noise(red_amp, red_gamma, fobs)[np.newaxis,:] # (1,F,) + noise = noise + red_noise # (P,F,) - # add single source noise - noise = noise[:,:,np.newaxis] - if ss_noise: - noise = noise + _Sh_ss_noise(hc_ss, fobs) # (P, F, R) + # add single source noise + noise = noise[:,:,np.newaxis] + if ss_noise: + noise = noise + _Sh_ss_noise(hc_ss, fobs) # (P, F, R) mu_1B = _mean1_Bstatistic(noise, Gamma, Sh_bg, Sh0_bg) @@ -872,7 +883,7 @@ def _antenna_pattern_functions(m_hat, n_hat, Omega_hat, pi_hat): ######################## Noise Spectral Density ######################## -def _Sh_rest_noise(hc_ss, hc_bg, freqs): +def _Sh_rest_noise(hc_ss, hc_bg, freqs, nexcl=0): """ Calculate the noise spectral density contribution from all but the current single source. Parameters @@ -883,6 +894,9 @@ def _Sh_rest_noise(hc_ss, hc_bg, freqs): Characteristic strain from all but loudest source at each frequency. freqs : (F,) 1Darray Frequency bin centers. + exclude_loudest : int + Number of loudest single sources to exclude from hc_rest noise, in addition + to the source in question. Returns ------- @@ -891,10 +905,15 @@ def _Sh_rest_noise(hc_ss, hc_bg, freqs): Follows Eq. (45) in Rosado et al. 2015. """ - hc2_louds = np.sum(hc_ss**2, axis=2) # (F,R) - # subtract the single source from rest of loud sources and the background, for each single source - hc2_rest = hc_bg[:,:,np.newaxis]**2 + hc2_louds[:,:,np.newaxis] - hc_ss**2 # (F,R,L) - Sh_rest = hc2_rest / freqs[:,np.newaxis,np.newaxis]**3 /(12 * np.pi**2) # (F,R,L) + + if nexcl>0: + Sh_rest = cyutils.Sh_rest(hc_ss, hc_bg, freqs, nexcl) + else: + hc2_louds = np.sum(hc_ss**2, axis=2) # (F,R) + # subtract the single source from rest of loud sources and the background, for each single source + hc2_rest = hc_bg[:,:,np.newaxis]**2 + hc2_louds[:,:,np.newaxis] - hc_ss**2 # (F,R,L) + Sh_rest = hc2_rest / freqs[:,np.newaxis,np.newaxis]**3 /(12 * np.pi**2) # (F,R,L) + return Sh_rest @@ -913,7 +932,7 @@ def _Sh_ss_noise(hc_ss, freqs): Returns ------- - Sh_ss : (F,R,L) NDarray of scalars + Sh_ss : (P,F,R) NDarray of scalars The noise in a single pulsar from other GW sources for detecting each single source. Follows Eq. (45) in Rosado et al. 2015. @@ -952,7 +971,8 @@ def _red_noise(red_amp, red_gamma, freqs, f_ref=1/YR): -def _total_noise(delta_t, sigmas, hc_ss, hc_bg, freqs, red_amp=None, red_gamma=None): +def _total_noise(delta_t, sigmas, hc_ss, hc_bg, freqs, red_amp=None, red_gamma=None, + nexcl=0): """ Calculate the noise spectral density of each pulsar, as it pertains to single source detections, i.e., including the background as a noise source. @@ -968,6 +988,8 @@ def _total_noise(delta_t, sigmas, hc_ss, hc_bg, freqs, red_amp=None, red_gamma=N Characteristic strain from all but loudest source at each frequency. freqs : (F,) 1Darray Frequency bin centers. + exclude_loudest : int + Number of loudest single sources to exclude from hc_rest noise Returns ------- @@ -978,7 +1000,7 @@ def _total_noise(delta_t, sigmas, hc_ss, hc_bg, freqs, red_amp=None, red_gamma=N """ noise = _white_noise(delta_t, sigmas) # (P,) - Sh_rest = _Sh_rest_noise(hc_ss, hc_bg, freqs) # (F,R,L,) + Sh_rest = _Sh_rest_noise(hc_ss, hc_bg, freqs, nexcl) # (F,R,L,) noise = noise[:,np.newaxis,np.newaxis,np.newaxis] + Sh_rest[np.newaxis,:,:,:] # (P,F,R,L) if (red_amp is not None) and (red_gamma is not None): red_noise = _red_noise(red_amp, red_gamma, freqs) # (F,) @@ -986,6 +1008,51 @@ def _total_noise(delta_t, sigmas, hc_ss, hc_bg, freqs, red_amp=None, red_gamma=N return noise +def psrs_spectra_gwbnoise(psrs, fobs, nreals, npsrs, divide_flag=False): + """ Get GWBSensitivityCurve noise and spectra for psrs + + Returns + ------- + spectra : + noise_gsc : (P,F,R) + """ + spectra = [] + for psr in psrs: + sp = hsen.Spectrum(psr, freqs=fobs) + sp.NcalInv + spectra.append(sp) + sc_bg = hsen.GWBSensitivityCurve(spectra).h_c + noise_gsc = sc_bg**2 / (12 *np.pi**2 *fobs**3) + noise_gsc = np.repeat(noise_gsc, npsrs*nreals).reshape(len(fobs), npsrs, nreals) # (F,P,R) + noise_gsc = np.swapaxes(noise_gsc, 0, 1) # (P,F,R) + + if divide_flag: noise_gsc *= npsrs*(npsrs-1) + + return spectra, noise_gsc + +def _dsc_noise(fobs, nreals, npsrs, nloudest, psrs=None, spectra=None, divide_flag=False): + """ Get DeterSensitivityCurve noise using either psrs or spectra + + Returns + ------- + noise_dsc : (P,F,R,L) NDarray + """ + + if spectra is None: + assert psrs is not None, 'Must provide spectra or psrs' + spectra = [] + for psr in psrs: + sp = hsen.Spectrum(psr, freqs=fobs) + sp.NcalInv + spectra.append(sp) + sc_ss = hsen.DeterSensitivityCurve(spectra).h_c + noise_dsc = sc_ss**2 / (12 *np.pi**2 *fobs**3) + noise_dsc = np.repeat(noise_dsc, npsrs*nreals*nloudest).reshape(len(fobs), npsrs, nreals, nloudest) # (F,P,R,L) + noise_dsc = np.swapaxes(noise_dsc, 0, 1) # (P,F,R,L) + + if divide_flag: noise_dsc *= npsrs*(npsrs-1) + return noise_dsc + ################### GW polarization, phase, amplitude ################### @@ -1093,7 +1160,7 @@ def _snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): """ - snr_ss = sam_cyutils.snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs) + snr_ss = cyutils.snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs) return snr_ss def _snr_ss_5dim(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): @@ -1128,24 +1195,16 @@ def _snr_ss_5dim(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): """ amp = amp[np.newaxis,:,:,np.newaxis,:] # (F,R,L) to (P,F,R,S,L) - # print('amp', amp.shape) a_pol, b_pol = _a_b_polarization(iotas) # (F,S,L) a_pol = a_pol[np.newaxis,:,np.newaxis,:,:] # (F,S,L) to (1,F,1,S,L) b_pol = b_pol[np.newaxis,:,np.newaxis,:,:] # (F,S,L) to (1,F,1,S,L) - # print('a_pol', a_pol.shape) - # print('b_pol', b_pol.shape) Phi_T = _gw_phase(dur, freqs, Phi_0) # (F,) - # print('Phi_T', Phi_T.shape) Phi_T = Phi_T[np.newaxis,:,np.newaxis,:,:] # (F,S,L) to (1,F,1,S,L) - # print('Phi_T', Phi_T.shape) Phi_0 = Phi_0[np.newaxis,:,np.newaxis,:,:] # (F,S,L) to (1,F,1,S,L) - # print('Phi_0', Phi_0.shape) - freqs = freqs[np.newaxis,:,np.newaxis,np.newaxis,np.newaxis] # (F,) to (1,F,1,1,1) - # print('freqs', freqs.shape) S_i = S_i[:,:,:,np.newaxis,:] # (P,F,R,L) to (P,F,R,1,L) F_iplus = F_iplus[:,:,np.newaxis,:,:] # (P,F,S,L) to (P,F,1,S,L) @@ -1347,14 +1406,14 @@ def _gamma_ssi_cython(rho, grid_path): # # interpolate for gamma in cython # rho_flat = rho[ff,rr].flatten() # rsort = np.argsort(rho_flat) - # gamma_flat = sam_cyutils.gamma_of_rho_interp(rho_flat, rsort, rho_interp_grid, gamma_interp_grid) + # gamma_flat = cyutils.gamma_of_rho_interp(rho_flat, rsort, rho_interp_grid, gamma_interp_grid) # gamma_ssi[ff,rr] = gamma_flat.reshape(rho[ff,rr].shape) for rr in range(len(rho[0])): # interpolate for gamma in cython rho_flat = rho[:,rr].flatten() rsort = np.argsort(rho_flat) - gamma_flat = sam_cyutils.gamma_of_rho_interp(rho_flat, rsort, rho_interp_grid, gamma_interp_grid) + gamma_flat = cyutils.gamma_of_rho_interp(rho_flat, rsort, rho_interp_grid, gamma_interp_grid) gamma_ssi[:,rr] = gamma_flat.reshape(rho[:,rr].shape) @@ -1439,6 +1498,8 @@ def detect_ss(thetas, phis, sigmas, fobs, hc_ss, hc_bg, SNR of each single source. """ + + warnings.warn("'detect_ss()' is deprecated. Use 'detect_ss_pta()' instead.", DeprecationWarning) dur = 1.0/fobs[0] cad = 1.0/(2*fobs[-1]) fobs_cents, fobs_edges = utils.pta_freqs(dur, num=len(fobs)) @@ -1485,7 +1546,8 @@ def detect_ss(thetas, phis, sigmas, fobs, hc_ss, hc_bg, return gamma_ss -def detect_ss_pta(pulsars, fobs, hc_ss, hc_bg, +def detect_ss_pta(pulsars, fobs, hc_ss, hc_bg, + custom_noise=None, nexcl_noise=0, theta_ss=None, phi_ss=None, Phi0_ss=None, iota_ss=None, psi_ss=None, nskies=25, Fe_bar = None, red_amp=None, red_gamma=None, alpha_0=0.001, Fe_bar_guess=15, ret_snr=False, print_nans=False, snr_cython=True, gamma_cython=True, grid_path=GAMMA_RHO_GRID_PATH): @@ -1504,6 +1566,11 @@ def detect_ss_pta(pulsars, fobs, hc_ss, hc_bg, hc_bg : (F,R) Characteristic strain of the background at each frequency, for R realizations. + custom_noise : (P,F,R,L) array or None + Custom noise if not None, otherwise noise is calcualted normally. + nexcl_noise : int + Number of loudest single sources to exclude from hc_rest noise, in addition + to the source in question. theta_ss : (F,S,L) NDarray Polar (latitudinal) angular position in the sky of each single source. Must be provided, to give the shape for sky realizations. @@ -1582,7 +1649,13 @@ def detect_ss_pta(pulsars, fobs, hc_ss, hc_bg, pi_hat) # (P,F,S,L) # noise spectral density - S_i = _total_noise(cad, sigmas, hc_ss, hc_bg, fobs, red_amp, red_gamma) + if custom_noise is not None: + if custom_noise.shape != (len(pulsars), nfreqs, nreals, nloudest): + err = f"{custom_noise.shape=}, must be shape (P,F,R,L)=({len(pulsars)}, {nfreqs}, {nreals}, {nloudest})" + raise ValueError(err) + S_i = custom_noise + else: + S_i = _total_noise(cad, sigmas, hc_ss, hc_bg, fobs, red_amp, red_gamma, nexcl=nexcl_noise) # amplitudw amp = _amplitude(hc_ss, fobs, dfobs) # (F,R,L) @@ -1775,8 +1848,8 @@ def detect_lib(hdf_name, output_dir, npsrs, sigma, nskies, thresh=DEF_THRESH, def detect_lib_clbrt_pta(hdf_name, output_dir, npsrs, nskies, thresh=DEF_THRESH, sigstart=1e-6, sigmin=1e-9, sigmax=1e-4, tol=0.01, maxbads=5, - plot=True, debug=False, grid_path=GAMMA_RHO_GRID_PATH, - snr_cython = True, save_ssi=False, ret_dict=False, ss_noise=True): + plot=True, debug=False, + save_ssi=False, ret_dict=False, ss_noise=False): """ Calculate detection statistics for an ss library output. Parameters @@ -1841,9 +1914,6 @@ def detect_lib_clbrt_pta(hdf_name, output_dir, npsrs, nskies, thresh=DEF_THRESH, fobs = ssfile['fobs'][:] dur = 1.0/fobs[0] cad = 1.0/(2*fobs[-1]) - # if dfobs is None: dfobs = ssfile['dfobs'][:] - # if dur is None: dur = ssfile['pta_dur'][0] - # if cad is None: cad = ssfile['pta_cad'][0] hc_ss = ssfile['hc_ss'][...] hc_bg = ssfile['hc_bg'][...] shape = hc_ss.shape @@ -1880,10 +1950,6 @@ def detect_lib_clbrt_pta(hdf_name, output_dir, npsrs, nskies, thresh=DEF_THRESH, snr_ss = np.zeros((nsamps, nfreqs, nreals, nskies, nloudest)) gamma_ssi = np.zeros((nsamps, nfreqs, nreals, nskies, nloudest)) - # # one time calculations - # Num = nfreqs * nloudest # number of single sources in a single strain realization (F*L) - # Fe_bar = _Fe_thresh(Num) # scalar - for nn in range(nsamps): if debug: print('\non sample nn=%d out of N=%d' % (nn,nsamps)) @@ -1902,7 +1968,7 @@ def detect_lib_clbrt_pta(hdf_name, output_dir, npsrs, nskies, thresh=DEF_THRESH, # use sigmin and sigmax from previous realization, # unless it's the first realization of the sample - psrs, _sigstart, _sigmin, _sigmax = calibrate_one_pta(hc_bg[nn,:,rr], fobs, npsrs, tol=tol, maxbads=maxbads, + psrs, _sigstart, _sigmin, _sigmax = calibrate_one_pta(hc_bg[nn,:,rr], hc_ss[nn,:,rr,:], fobs, npsrs, tol=tol, maxbads=maxbads, sigstart=_sigstart, sigmin=_sigmin, sigmax=_sigmax, debug=debug, ret_sig=True, ss_noise=ss_noise) _sigmin /= 2 _sigmax *= 2 @@ -1926,15 +1992,7 @@ def detect_lib_clbrt_pta(hdf_name, output_dir, npsrs, nskies, thresh=DEF_THRESH, now = datetime.now() print(f"Sample {nn} took {now-samp_dur} s") samp_dur = now - # dp_bg[nn,:], snr_bg[nn,...] = detect_bg_pta(psrs, fobs, hc_bg[nn], ret_snr=True) - # vals_ss = detect_ss_pta(psrs, fobs, hc_ss[nn], hc_bg[nn], - # ret_snr=True, gamma_cython=True, snr_cython=snr_cython, - # theta_ss=theta_ss, phi_ss=phi_ss, Phi0_ss=Phi0_ss, - # iota_ss=iota_ss, psi_ss=psi_ss, grid_path=grid_path) - # dp_ss[nn,:,:] = vals_ss[0] - # if save_ssi: - # snr_ss[nn] = vals_ss[1] - # gamma_ssi[nn] = vals_ss[2] + df_ss[nn], df_bg[nn] = detfrac_of_reals(dp_ss[nn], dp_bg[nn], thresh) if save_ssi: @@ -2025,7 +2083,7 @@ def detfrac_of_reals(dp_ss, dp_bg, thresh=DEF_THRESH): def expval_of_ss(gamma_ssi,): - """ Calculate the expected number of single source detections, across all realization + """ Calculate the expected number of single source detections, across all frequencies. Parameters ---------- @@ -2035,7 +2093,7 @@ def expval_of_ss(gamma_ssi,): Returns ------- ev_ss : (R,S) - Expected number of single source detection (dp_ss>thresh) for each strain and sky realizations. + Expected (fractional) number of single source detection for each strain and sky realizations. """ # print(f"{gamma_ssi.shape=}, {[*gamma_ssi.shape]}") @@ -2044,6 +2102,26 @@ def expval_of_ss(gamma_ssi,): return ev_ss # df_bg[nn] = np.sum(dp_bg[nn]>thresh)/(nreals) +def count_n_ss(gamma_ssi): + """ Calculate the number of random single source detections. + + Parameters + ---------- + gamma_ssi : (F,R,S,L) NDarray + Detection probability of each single source + + Returns + ------- + nn_ss : (R,S) + Number of random single source detections for each strain and sky realization. + + """ + + randoms = np.random.uniform(0,1, size=gamma_ssi.size).reshape(gamma_ssi.shape) + nn_ss = np.sum((gamma_ssi>randoms), axis=(0,3)) + return nn_ss + + ############################# Plot Library ############################# @@ -2277,24 +2355,88 @@ def detect_pspace_model(fobs_cents, hc_ss, hc_bg, return dsdata -def detect_pspace_model_clbrt_pta(fobs_cents, hc_ss, hc_bg, npsrs, nskies, - sigstart=1e-6, sigmin=1e-9, sigmax=1e-4, tol=0.01, maxbads=5, - thresh=DEF_THRESH, debug=False, save_snr_ss=False, save_gamma_ssi=True, - red_amp=None, red_gamma=None, red2white=None, ss_noise=True): +def detect_pspace_model_psrs(fobs_cents, hc_ss, hc_bg, psrs, nskies, hc_bg_noise=None, + thresh=DEF_THRESH, debug=False, nexcl_noise=0): + + nfreqs, nreals, nloudest = [*hc_ss.shape] + dur = 1/fobs_cents[0] + cad = 1.0 / (2 * fobs_cents[-1]) + # fobs_cents, fobs_edges = holo.utils.pta_freqs(dur) + # dfobs = np.diff(fobs_edges) + + if hc_bg_noise is None: + hc_bg_noise=hc_bg + elif debug: + print("Using different hc_bg_noise for SS detstats.") + + + # Build ss skies + if debug: print('Building ss skies.') + theta_ss, phi_ss, Phi0_ss, iota_ss, psi_ss = _build_skies(nfreqs, nskies, nloudest) + + + # Calculate DPs, SNRs, and DFs + if debug: print('Calculating SS and BG detection statistics.') + dp_bg, snr_bg = detect_bg_pta(psrs, fobs_cents, hc_bg, ret_snr=True) + # print(f"{np.mean(dp_bg)=}") + vals_ss = detect_ss_pta( + psrs, fobs_cents, hc_ss, hc_bg_noise, + gamma_cython=True, snr_cython=True, ret_snr=True, + theta_ss=theta_ss, phi_ss=phi_ss, Phi0_ss=Phi0_ss, iota_ss=iota_ss, psi_ss=psi_ss, + nexcl_noise=nexcl_noise + ) + dp_ss, snr_ss, gamma_ssi = vals_ss[0], vals_ss[1], vals_ss[2] + # print(f"{np.mean(dp_ss)=}") + df_ss = np.sum(dp_ss>thresh)/(nreals*nskies) + df_bg = np.sum(dp_bg>thresh)/(nreals) + ev_ss = expval_of_ss(gamma_ssi) + # print(f"{np.mean(ev_ss)=}") + + dsdata = { + 'dp_ss':dp_ss, 'snr_ss':snr_ss, 'gamma_ssi':gamma_ssi, + 'dp_bg':dp_bg, 'snr_bg':snr_bg, + 'df_ss':df_ss, 'df_bg':df_bg, 'ev_ss':ev_ss, + } + + return dsdata + +def detect_pspace_model_clbrt_pta( + fobs_cents, hc_ss, hc_bg, npsrs, nskies, + hc_bg_noise=None, + sigstart=1e-6, sigmin=1e-9, sigmax=1e-4, tol=0.01, maxbads=5, + thresh=DEF_THRESH, debug=False, save_snr_ss=False, save_gamma_ssi=True, + red_amp=None, red_gamma=None, red2white=None, ss_noise=False, dsc_flag=False, nexcl_noise=0): """ Detect pspace model using individual sigma calibration for each realization Parameters ---------- - + fobs_cents : 1Darray + GW frequencies in Hz + hc_ss : (F,R,L) NDarray + hc_bg : (F,R) NDarray + npsrs : int + nskies : int + hc_bg_noise " (F,R) NDarray or None + the background to use as single source noise, when normal hc_bg has extra SS added to it. red2white : scalar or None Fixed ratio between red and white noise amplitude, if not None. Otherwise, red noise stays fixed - + nexcl_noise : int + Number of loudest single sources to exclude from hc_rest noise, in addition + to the source in question. """ dur = 1.0/fobs_cents[0] cad = 1.0/(2*fobs_cents[-1]) nfreqs, nreals, nloudest = [*hc_ss.shape] + + if hc_bg_noise is None: + hc_bg_noise=hc_bg + elif debug: + print("Using different hc_bg_noise for SS detstats.") + if np.any(hc_bg_noise>hc_bg): + err = f"hc_bg_noise excluding SS is somehow larger than hc_bg with SS added in!!" + raise ValueError(err) # form arrays for individual realization detstats # set all to nan, only to be replaced if successful pta is found @@ -2303,6 +2445,7 @@ def detect_pspace_model_clbrt_pta(fobs_cents, hc_ss, hc_bg, npsrs, nskies, snr_ss = np.ones((nfreqs, nreals, nskies, nloudest)) * np.nan snr_bg = np.ones((nreals)) * np.nan gamma_ssi = np.ones((nfreqs, nreals, nskies, nloudest)) * np.nan + sigmas = np.ones((nreals)) * np.nan # for each realization, @@ -2320,28 +2463,161 @@ def detect_pspace_model_clbrt_pta(fobs_cents, hc_ss, hc_bg, npsrs, nskies, print(f"{rr=}, {(now-real_dur)/100} s per realization, {_sigmin=:.2e}, {_sigmax=:.2e}, {_sigstart=:.2e}") real_dur = now + # get calibrated psrs psrs, red_amp, _sigstart, _sigmin, _sigmax = calibrate_one_pta(hc_bg[:,rr], hc_ss[:,rr,:], fobs_cents, npsrs, tol=tol, maxbads=maxbads, sigstart=_sigstart, sigmin=_sigmin, sigmax=_sigmax, debug=debug, ret_sig=True, red_amp=red_amp, red_gamma=red_gamma, red2white=red2white, ss_noise=ss_noise) _sigmin /= 2 _sigmax *= 2 + 2e-20 # >1e-20 to make sure it doesnt immediately fail the 0 check + sigmas[rr] = _sigstart if psrs is None: failed_psrs += 1 continue # leave values as nan, if no successful PTA was found - # print(f"before calculation: {utils.stats(psrs[0].toaerrs)=}, \n{utils.stats(hc_bg[rr])=},\ - # {utils.stats(fobs_cents)=}") + + # use those psrs to calculate realization detstats _dp_bg, _snr_bg = detect_bg_pta(psrs, fobs_cents, hc_bg[:,rr:rr+1], hc_ss[:,rr:rr+1,:], ret_snr=True, red_amp=red_amp, red_gamma=red_gamma) - # print(f"{utils.stats(psrs[0].toaerrs)=}, {utils.stats(hc_bg[rr])=},\ - # {_dp_bg=},") - # _dp_bg, = detect_bg_pta(psrs, fobs_cents, hc_bg=hc_bg[:,rr:rr+1], red_amp=red_amp, red_gamma=red_gamma) #, ret_snr=True) - # print(f"test2: {_dp_bg=}") + dp_bg[rr], snr_bg[rr] = _dp_bg.squeeze(), _snr_bg.squeeze() + + + # calculate SS noise from DeterSensitivityCurve and S_h,rest + if dsc_flag: + spectra = [] + for psr in psrs: + sp = hsen.Spectrum(psr, freqs=fobs_cents) + sp.NcalInv + spectra.append(sp) + sc_hc = hsen.DeterSensitivityCurve(spectra).h_c + noise_dsc = sc_hc**2 / (12 * np.pi**2 * fobs_cents**3) + noise_dsc = _dsc_noise(fobs_cents, nreals, npsrs, nloudest, psrs, spectra) # (P,F,R,L) + # np.repeat(noise_dsc, npsrs*1*nloudest).reshape(nfreqs, npsrs, 1, nloudest) # (F,P,R,L) + # noise_dsc = np.swapaxes(noise_dsc, 0, 1) # (P,F,R,L) + noise_rest = _Sh_rest_noise(hc_ss[:,rr:rr+1,:], hc_bg[:,rr:rr+1], fobs_cents) # (F,R,L) + noise_ss = noise_dsc + noise_rest[np.newaxis,:,:,:] # (P,F,R,L) + else: + noise_ss = None + + # calculate realization SS detstats _dp_ss, _snr_ss, _gamma_ssi = detect_ss_pta( - psrs, fobs_cents, hc_ss[:,rr:rr+1], hc_bg[:,rr:rr+1], nskies=nskies, ret_snr=True, red_amp=red_amp, red_gamma=red_gamma) - # if debug: print(f"{_dp_ss.shape=}, {_snr_ss.shape=}, {_gamma_ssi.shape=}") + psrs, fobs_cents, hc_ss[:,rr:rr+1], hc_bg_noise[:,rr:rr+1], custom_noise=noise_ss, + nskies=nskies, ret_snr=True, red_amp=red_amp, red_gamma=red_gamma, nexcl_noise=nexcl_noise) + + dp_ss[rr] = _dp_ss.reshape(nskies) # from R=1,S to S + snr_ss[:,rr] = _snr_ss.reshape(nfreqs, nskies, nloudest) # from F,R=1,S,L to F,S,L + gamma_ssi[:,rr,:,:] = _gamma_ssi.reshape(nfreqs, nskies, nloudest) # from F,R=1,S,L to F,S,L + + ev_ss = expval_of_ss(gamma_ssi) + df_ss, df_bg = detfrac_of_reals(dp_ss, dp_bg) + _dsdat = { + 'dp_ss':dp_ss, 'snr_ss':snr_ss, 'gamma_ssi':gamma_ssi, + 'dp_bg':dp_bg, 'snr_bg':snr_bg, + 'df_ss':df_ss, 'df_bg':df_bg, 'ev_ss':ev_ss, + 'sigmas':sigmas, + } + if save_gamma_ssi: + _dsdat.update(gamma_ssi=gamma_ssi) + if save_snr_ss: + _dsdat.update(snr_ss=snr_ss) + if debug: + print(f"Model took {datetime.now() - mod_start} s, {failed_psrs}/{nreals} realizations failed.") + return _dsdat + + +def detect_pspace_model_clbrt_pta_gsc( + fobs_cents, hc_ss, hc_bg, npsrs, nskies, hc_bg_noise=None, + sigstart=1e-6, sigmin=1e-9, sigmax=1e-4, tol=0.01, maxbads=5, + thresh=DEF_THRESH, debug=False, save_snr_ss=False, save_gamma_ssi=True, + red_amp=None, red_gamma=None, red2white=None, ss_noise=False, + divide_flag=False, dsc_flag=False, nexcl_noise=0): + """ Detect pspace model using individual sigma calibration for each realization + and sensitivity curve noise for both BG calibration and SS detstats. + + Parameters + ---------- + fobs_cents : 1Darray + GW frequencies in Hz + hc_ss : (F,R,L) NDarray + hc_bg : (F,R) NDarray + npsrs : int + nskies : int + hc_bg_noise " (F,R) NDarray or None + the background to use as single source noise, when normal hc_bg has extra SS added to it. + red2white : scalar or None + Fixed ratio between red and white noise amplitude, if not None. + Otherwise, red noise stays fixed + nexcl_noise : int + Number of loudest single sources to exclude from hc_rest noise, in addition + to the source in question. + """ + dur = 1.0/fobs_cents[0] + cad = 1.0/(2*fobs_cents[-1]) + if hc_bg_noise is None: + hc_bg_noise=hc_bg + + nfreqs, nreals, nloudest = [*hc_ss.shape] + + # form arrays for individual realization detstats + # set all to nan, only to be replaced if successful pta is found + dp_ss = np.ones((nreals, nskies)) * np.nan + dp_bg = np.ones(nreals) * np.nan + snr_ss = np.ones((nfreqs, nreals, nskies, nloudest)) * np.nan + snr_bg = np.ones((nreals)) * np.nan + gamma_ssi = np.ones((nfreqs, nreals, nskies, nloudest)) * np.nan + + + # for each realization, + # use sigmin and sigmax from previous realization, + # unless it's the first realization of the sample + _sigstart, _sigmin, _sigmax = sigstart, sigmin, sigmax + if debug: + mod_start = datetime.now() + real_dur = datetime.now() + failed_psrs=0 + for rr in range(nreals): + if debug: + now = datetime.now() + if (rr%100==99): + print(f"{rr=}, {(now-real_dur)/100} s per realization, {_sigmin=:.2e}, {_sigmax=:.2e}, {_sigstart=:.2e}") + real_dur = now + + + # get calibrated psrs + psrs, red_amp, _sigstart, _sigmin, _sigmax, spectra, noise_gsc = calibrate_one_pta_gsc( + hc_bg[:,rr], fobs_cents, npsrs, tol=tol, maxbads=maxbads, + sigstart=_sigstart, sigmin=_sigmin, sigmax=_sigmax, debug=debug, ret_sig=True, + red_amp=red_amp, red_gamma=red_gamma, red2white=red2white, ss_noise=ss_noise, divide_flag=divide_flag) + _sigmin /= 2 + _sigmax *= 2 + 2e-20 # >1e-20 to make sure it doesnt immediately fail the 0 check + + if psrs is None: + failed_psrs += 1 + continue # leave values as nan, if no successful PTA was found + + + # use those psrs to calculate realization detstats + _dp_bg, _snr_bg = detect_bg_pta(psrs, fobs_cents, hc_bg[:,rr:rr+1], hc_ss[:,rr:rr+1,:], custom_noise=noise_gsc, + ret_snr=True, red_amp=red_amp, red_gamma=red_gamma) + + dp_bg[rr], snr_bg[rr] = _dp_bg.squeeze(), _snr_bg.squeeze() + + + # calculate SS noise from DeterSensitivityCurve and S_h,rest + if dsc_flag: # calculate the DSC noise separately + noise_dsc = _dsc_noise(fobs_cents, spectra=spectra, nreals=1, npsrs=npsrs, nloudest=nloudest, divide_flag=divide_flag) + else: # use the GSC as the DSC noise + noise_dsc = np.repeat(noise_gsc, nloudest).reshape(npsrs, nfreqs, 1, nloudest) # (P,F,R,L) + noise_rest = _Sh_rest_noise(hc_ss[:,rr:rr+1,:], hc_bg[:,rr:rr+1], fobs_cents) # (F,R,L) + noise_ss = noise_dsc + noise_rest[np.newaxis,:,:,:] # (P,F,R,L) + + + # calculate realizatoin SS detstats + _dp_ss, _snr_ss, _gamma_ssi = detect_ss_pta( + psrs, fobs_cents, hc_ss[:,rr:rr+1], hc_bg_noise[:,rr:rr+1], custom_noise=noise_ss, + nskies=nskies, ret_snr=True, red_amp=red_amp, red_gamma=red_gamma, nexcl_noise=nexcl_noise) + dp_ss[rr], snr_ss[:,rr], gamma_ssi[:,rr] = _dp_ss.squeeze(), _snr_ss.squeeze(), _gamma_ssi.squeeze() ev_ss = expval_of_ss(gamma_ssi) @@ -2362,7 +2638,7 @@ def detect_pspace_model_clbrt_pta(fobs_cents, hc_ss, hc_bg, npsrs, nskies, def detect_pspace_model_clbrt_ramp(fobs_cents, hc_ss, hc_bg, npsrs, nskies, sigma, rampstart=1e-16, rampmin=1e-20, rampmax=1e-13, tol=0.01, maxbads=5, thresh=DEF_THRESH, debug=False, save_snr_ss=False, save_gamma_ssi=True, - red_amp=None, red_gamma=None, ss_noise=True): + red_amp=None, red_gamma=None, ss_noise=False): """ Detect pspace model using individual red noise amplitude calibration for each realization NOTE: Not supported, not updated for including single sources as noise for BG. @@ -2572,10 +2848,16 @@ def calibrate_all_sigma(hc_bg, fobs, npsrs, maxtrials, ) return rsigmas, avg_dps, std_dps + +def _red_amp_from_white_noise(cad, sigma, red2white, fref=1/YR): + red_amp = np.sqrt(12 * np.pi**2 * fref**3 * + red2white * _white_noise(cad, sigma)) + return red_amp + def calibrate_one_pta(hc_bg, hc_ss, fobs, npsrs, sigstart=1e-6, sigmin=1e-9, sigmax=1e-4, debug=False, maxbads=20, tol=0.03, phis=None, thetas=None, ret_sig = False, red_amp=None, red_gamma=None, red2white=None, - ss_noise=True): + ss_noise=False): """ Calibrate the specific PTA for a given realization, and return that PTA Parameters @@ -2611,7 +2893,7 @@ def calibrate_one_pta(hc_bg, hc_ss, fobs, npsrs, if thetas is None: thetas = np.random.uniform(np.pi/2, np.pi/2, size = npsrs) sigma = sigstart if red2white is not None: - red_amp = _white_noise(cad, sigma) * red2white + red_amp = _red_amp_from_white_noise(cad, sigma, red2white) psrs = hsim.sim_pta(timespan=dur/YR, cad=1/(cad/YR), sigma=sigma, phi=phis, theta=thetas) @@ -2625,7 +2907,7 @@ def calibrate_one_pta(hc_bg, hc_ss, fobs, npsrs, while np.abs(dp_bg-0.50)>tol: sigma = np.mean([sigmin, sigmax]) # a weighted average would be better if red2white is not None: - red_amp = sigma * red2white + red_amp = _red_amp_from_white_noise(cad, sigma, red2white) psrs = hsim.sim_pta(timespan=dur/YR, cad=1/(cad/YR), sigma=sigma, phi=phis, theta=thetas) dp_bg = detect_bg_pta(psrs, fobs, hc_bg=hc_bg[:,np.newaxis], hc_ss=hc_ss[:,np.newaxis,:], @@ -2676,11 +2958,126 @@ def calibrate_one_pta(hc_bg, hc_ss, fobs, npsrs, # {utils.stats(fobs)=}, {dp_bg=}") if ret_sig: return psrs, red_amp, sigma, sigmin, sigmax + return psrs + + +def calibrate_one_pta_gsc(hc_bg, fobs, npsrs, + sigstart=1e-6, sigmin=1e-9, sigmax=1e-4, debug=False, maxbads=20, tol=0.03, + phis=None, thetas=None, ret_sig = False, red_amp=None, red_gamma=None, red2white=None, + ss_noise=False, divide_flag=False,): + """ Calibrate the specific PTA for a given realization, and return that PTA + + Parameters + ---------- + hc_bg : (F,) 1Darray + The background characteristic strain for one realization. + hc_ss : (F,L) NDarray + The SS characteristic strains for one realization + fobs : (F,) 1Darray + Observed GW frequencies. + npsrs : integer + Number of pulsars. + divide_flag : Bool + Whether or not to divide the GSC noise among the pulsars + + Returns + ------- + psrs : hasasia.sim.pta object + Calibrated PTA. + sigmin : float + minimum of the final sigma range used, returned only if ret_sig=True + sigmax : float, returned only if ret_sig=True + maximum of the final sigma range used + sigma : float + final sigma, returned only if ret_sig=True + + """ + + # get duration and cadence from fobs + dur = 1.0/fobs[0] + cad = 1.0/(2.0*fobs[-1]) + + # randomize pulsar positions + if phis is None: phis = np.random.uniform(0, 2*np.pi, size = npsrs) + if thetas is None: thetas = np.random.uniform(np.pi/2, np.pi/2, size = npsrs) + sigma = sigstart + if red2white is not None: + red_amp = _white_noise(cad, sigma) * red2white + + psrs = hsim.sim_pta(timespan=dur/YR, cad=1/(cad/YR), sigma=sigma, + phi=phis, theta=thetas) + + # get sensitivity curve + spectra, noise_gsc = psrs_spectra_gwbnoise(psrs, fobs, nreals=1, npsrs=npsrs, divide_flag=divide_flag) + + dp_bg = detect_bg_pta(psrs, fobs, hc_bg=hc_bg[:,np.newaxis], custom_noise=noise_gsc, + red_amp=red_amp, red_gamma=red_gamma, ss_noise=ss_noise)[0] + + nclose=0 # number of attempts close to 0.5, could be stuck close + nfar=0 # number of attempts far from 0.5, could be stuck far + + # calibrate sigma + while np.abs(dp_bg-0.50)>tol: + sigma = np.mean([sigmin, sigmax]) # a weighted average would be better + if red2white is not None: + red_amp = sigma * red2white + psrs = hsim.sim_pta(timespan=dur/YR, cad=1/(cad/YR), sigma=sigma, + phi=phis, theta=thetas) + spectra, noise_gsc = psrs_spectra_gwbnoise(psrs, fobs, nreals=1, npsrs=npsrs) + dp_bg = detect_bg_pta(psrs, fobs, hc_bg=hc_bg[:,np.newaxis], custom_noise=noise_gsc, + red_amp=red_amp, red_gamma=red_gamma, ss_noise=ss_noise)[0] + + # if debug: print(f"{dp_bg=}") + if (dp_bg < (0.5-tol)) or (dp_bg > (0.5+tol)): + nfar +=1 + + # check if we need to expand the range + if (nfar>5*maxbads # if we've had many bad guesses + or (sigmin/sigmax > 0.99 # or our range is small and we are far from the goal + and (dp_bg<0.4 or dp_bg>0.6))): + + # then we must expand the range + if debug: print(f"STUCK! {nfar=}, {dp_bg=}, {sigmin=:e}, {sigmax=:e}") + if dp_bg < 0.5-tol: # stuck way too low, allow much lower sigmin to raise DP + sigmin = sigmin/3 + if dp_bg > 0.5+tol: # stuck way too high, allow much higher sigmax to lower DP + sigmax = sigmax*3 + + # reset count for far guesses + nfar = 0 + + # check how we should narrow our range + if dp_bg<0.5-tol: # dp too low, lower sigma + sigmax = sigma + elif dp_bg>0.5+tol: # dp too high, raise sigma + sigmin = sigma + else: + nclose += 1 # check how many attempts between 0.49 and 0.51 fail + + # check if we are stuck near the goal value with a bad range + if nclose>maxbads: # if many fail, we're stuck; expand sampling range + if debug: print(f"{nclose=}, {dp_bg=}, {sigmin=:e}, {sigmax=:e}") + sigmin = sigmin/3 + sigmax = sigmax*3 + nclose=0 + + # check if goal DP is just impossible + if sigmax<1e-20: + psrs=None + if debug: print(f"FAILED! DP_BG=0.5 impossible with {red_amp=}, {red_gamma=}") + break + # print(f"test1: {dp_bg=}") + # print(f"test1: {sigma=}") + # print(f"in calibration: {utils.stats(psrs[0].toaerrs)=}, \n{utils.stats(hc_bg)=},\ + # {utils.stats(fobs)=}, {dp_bg=}") + if ret_sig: + return psrs, red_amp, sigma, sigmin, sigmax, spectra, noise_gsc return psrs, red_amp + def calibrate_one_ramp(hc_bg, hc_ss, fobs, psrs, rampstart=1e-6, rampmin=1e-9, rampmax=1e-4, debug=False, maxbads=20, tol=0.03, - phis=None, thetas=None, rgam=-1.5, ss_noise=True): + phis=None, thetas=None, rgam=-1.5, ss_noise=False): """ Calibrate the red noise amplitude, for a given realization, and return that PTA Parameters @@ -2768,4 +3165,661 @@ def calibrate_one_ramp(hc_bg, hc_ss, fobs, psrs, # print(f"test1: {sigma=}") # print(f"in calibration: {utils.stats(psrs[0].toaerrs)=}, \n{utils.stats(hc_bg)=},\ # {utils.stats(fobs)=}, {dp_bg=}") - return ramp, rampmin, rampmax \ No newline at end of file + return ramp, rampmin, rampmax + + + +######################################################################## +########################## Average Frequency ########################### +######################################################################## + +def weighted_mean_variance(data, weights, debug=False,): + """ Calculate the weighted average frequency and variance + + Parameters + ---------- + data: NDarray + Data + weights: NDarray + Weights + + Returns + ------- + mean : float + Weighted mean + var2 : float + Weighted variance (std^2) + """ + mean = np.sum(weights * data) / np.sum(weights) + if debug: print(f"{mean=}") + var2 = np.sum(weights * (data - mean)**2) + nn = data.size + var2 /= (nn-1)/nn * np.sum(weights) + if debug: print(f"{var2=}") + return mean, var2 + + + + + + + + +######################################################################## +###################### Functions Using V21 Models ###################### +######################################################################## + +def get_data( + target, dets=True, + nvars=21, nreals=500, nskies=100, shape=None, # keep as defaults + nloudest = 10, bgl = 10, cv=None, ssn_flag=False, + red_gamma = None, red2white=None, + gsc_flag=False, dsc_flag=False, divide_flag=False, nexcl=0, + gw_only=False, + var_hard_time=None, npsrs=40, +): + if gw_only: + path = '/Users/emigardiner/GWs/holodeck/output/anatomy_7GW' + else: + path = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz' + + output_path = path+f'/{target}_v{nvars}_r{nreals}_shape{str(shape)}' + + if var_hard_time is not None: + output_path+=f"_vtau{var_hard_time}" + + data_file = output_path +f'/data_params' + dets_file = output_path + f'/detstats_s{nskies}' + + + if nloudest != 10: # if using nloudest that isn't the default 10 + dets_file += f"_l{nloudest}" + data_file += f"_l{nloudest}" + if bgl != nloudest: + dets_file += f"_bgl{bgl}" # only change nloudest subtracted from bg, not single sources loudest + if cv is not None: + dets_file += f"_cv{cv}" # if using one variation to calibrate + if ssn_flag: + dets_file += '_ssn' # if using single sources as noise + + if gsc_flag: # if using GSC as noise + dets_file += '_gsc' + if dsc_flag is False: # if using GSC as noise + dets_file += 'both' + if divide_flag: + dets_file += '_divide' + else: + dets_file += '_nodiv' + if dsc_flag: + dets_file += '_dsc' + + if nexcl>0: + dets_file += f'_nexcl{nexcl}' + + if npsrs != 40: + dets_file += f'_p{npsrs}' + + if red2white is not None and red_gamma is not None: # if using red noise with fixed red_gamma + dets_file += f'_r2w{red2white:.1e}_rg{red_gamma:.1f}' + else: + dets_file += f'_white' + + dets_file += '.npz' + data_file += '.npz' + + print(f'{data_file=}') + print(f'{dets_file=}') + + if os.path.exists(data_file) is False: + err = f"load data file '{data_file}' does not exist, you need to construct it." + raise Exception(err) + if os.path.exists(dets_file) is False and dets is True: + err = f"load dets file '{dets_file}' does not exist, you need to construct it." + raise Exception(err) + file = np.load(data_file, allow_pickle=True) + data = file['data'] + params = file['params'] + file.close() + if dets is False: + return data, params + print(f"got data from {data_file}") + file = np.load(dets_file, allow_pickle=True) + print(f"got detstats from {dets_file}") + # print(file.files) + dsdat = file['dsdat'] + file.close() + + return data, params, dsdat + +def append_filename(filename='', + gw_only=False, red_gamma = None, red2white=None, + nloudest = 10, bgl = 10, cv=None, + gsc_flag=False, dsc_flag=False, divide_flag=False, nexcl=0, + var_hard_time=None, + + ): + + if cv is not None: + filename += f"_cv{cv}" + + if nloudest != 10: + filename += f"_l{nloudest}" + if bgl != nloudest: + filename += f"_bgl{bgl}" + + if var_hard_time is not None: + filename += f"_vtau{var_hard_time}" + + if red2white is not None and red_gamma is not None: + filename += f"_r2w{red2white:.1e}_rg{red_gamma:.1f}" + + if gsc_flag: + filename = filename + '_gsc' + if dsc_flag: filename += 'both' + if divide_flag: + filename += '_divide' + else: + filename += '_nodiv' + elif dsc_flag: filename += '_dsc' + + if nexcl>0: + filename += f'_nexcl{nexcl}' + + if gw_only: + filename = filename+'_gw' + + + return filename + +def build_ratio_arrays( + target, nreals=500, nskies=100, + gw_only=False, red2white=None, red_gamma=None, + nloudest=10, bgl=1, + gsc_flag=False, dsc_flag=False, divide_flag=False, nexcl=0, + var_hard_time=None, npsrs=40, + figpath = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/ratio',): + + data, params, dsdat = get_data(target, + gw_only=gw_only, red2white=red2white, red_gamma=red_gamma, + nloudest=nloudest, bgl=bgl, nexcl=nexcl, + gsc_flag=gsc_flag, dsc_flag=dsc_flag, divide_flag=divide_flag, + var_hard_time=var_hard_time, npsrs=npsrs) + xx=[] + yy=[] + for pp, par in enumerate(params): + xx.append(params[pp][target]) + dp_bg = np.repeat(dsdat[pp]['dp_bg'], nskies).reshape(nreals, nskies) + ev_ss = dsdat[pp]['ev_ss'] + yy.append(ev_ss/dp_bg) + + filename = figpath+f'/ratio_arrays_{target}' + filename = append_filename( + filename, + gw_only=gw_only, red_gamma=red_gamma, red2white=red2white, + nloudest=nloudest, bgl=bgl, cv=None, + gsc_flag=gsc_flag, dsc_flag=dsc_flag, divide_flag=divide_flag, nexcl=nexcl, + var_hard_time=var_hard_time) + filename += '.npz' + print(f'saving to {filename}') + np.savez(filename, xx_params = xx, yy_ratio = yy,) + +def get_ratio_arrays( + target, nreals=500, nskies=100, + gw_only=False, red2white=None, red_gamma=None, + nloudest=10, bgl=1, + gsc_flag=False, dsc_flag=False, divide_flag=False, nexcl=0, + var_hard_time=None, + figpath = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/ratio',): + + filename = figpath+f'/ratio_arrays_{target}' + filename = append_filename( + filename, + gw_only=gw_only, red_gamma=red_gamma, red2white=red2white, + nloudest=nloudest, bgl=bgl, cv=None, + gsc_flag=gsc_flag, dsc_flag=dsc_flag, divide_flag=divide_flag, nexcl=nexcl, + var_hard_time=var_hard_time) + filename += '.npz' + + file = np.load(filename) + xx_params = file['xx_params'] + yy_ratio = file['yy_ratio'] + file.close() + return xx_params, yy_ratio + +def build_noise_arrays( + target, nreals=500, nskies=100, + gw_only=False, red2white=None, red_gamma=None, + nloudest=10, bgl=1, save_temp=True, + gsc_flag=False, dsc_flag=False, divide_flag=False, nexcl=0, + var_hard_time=None, + figpath = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/noise',): + + data, params, dsdat = get_data(target, + gw_only=gw_only, red2white=red2white, red_gamma=red_gamma, + nloudest=nloudest, bgl=bgl, nexcl=nexcl, + gsc_flag=gsc_flag, dsc_flag=dsc_flag, divide_flag=divide_flag, + var_hard_time=var_hard_time) + fobs_cents = data[0]['fobs_cents'] + nfreqs=len(fobs_cents) + cad = 1.0/(2.0*fobs_cents[-1]) + + sigmas = [] + hc_ss = [] + hc_bg = [] + dp_max = [] + dp_2nd = [] + count_cws_50 = [] # number of single sources with DP>0.5 in any realization + count_cws_01 = [] # number of single sources with DP>0.5 in any realization + for ii, dat in enumerate(data): + sigmas.append(dsdat[ii]['sigmas']) # R, + hc_ss.append(dat['hc_ss']) + hc_bg.append(dat['hc_bg']) + + dp_ssi = dsdat[ii]['gamma_ssi'] # F,R,S,L + count_cws_01.append(np.sum(dp_ssi>0.01, axis=(0,3))) # from F,R,S,L to R,S + count_cws_50.append(np.sum(dp_ssi>0.50, axis=(0,3))) + + dp_ssi = np.swapaxes(dp_ssi, 1,3).reshape(nfreqs*nloudest, nreals*nskies) # F*L, S*R + argmax = np.argmax(dp_ssi, axis=0) # S*R + reals = np.arange(nreals*nskies) # S*R + _dp_max = dp_ssi[argmax, reals] # S*R + + dp_ssi[argmax, reals] = 0 + argmax = np.argmax(dp_ssi, axis=0) # S*R + _dp_2nd = dp_ssi[argmax, reals] # S*R + + dp_max.append(_dp_max) + dp_2nd.append(_dp_2nd) + + sigmas = np.array(sigmas) # V, R + hc_ss = np.array(hc_ss) # (V,F,R,L) + hc_bg = np.array(hc_bg) # (V,F,R) + count_cws_50 = np.array(count_cws_50) # V,R,S + count_cws_01 = np.array(count_cws_01) # V,R,S + dp_max=np.array(dp_max) # V,R,S, + dp_2nd=np.array(dp_2nd) # V,R,S, + + temp_name = '/Users/emigardiner/GWs/holodeck/output/temp' + temp_name += f'/noise_arrays_{target}.npz' + + white_noise = _white_noise(cad, sigmas) # V,R, array + if red2white is not None: + red_amp = _red_amp_from_white_noise(cad, sigmas, red2white) # V,R, array + red_noise = _red_noise(red_amp[np.newaxis,:,:], # (V,1,R,) + red_gamma, + fobs_cents[np.newaxis,:,np.newaxis] # (1,F,1) + ) + np.savez(temp_name, white_noise=white_noise, red_noise=red_noise, + count_cws_50=count_cws_50, count_cws_01=count_cws_01, + hc_ss=hc_ss, hc_bg=hc_bg) + return white_noise, red_noise, count_cws_50, count_cws_01, hc_ss, hc_bg, dp_max, dp_2nd + np.savez(temp_name, white_noise=white_noise, + count_cws_50=count_cws_50, count_cws_01=count_cws_01, + hc_ss=hc_ss, hc_bg=hc_bg, dp_max=dp_max, dp_2nd=dp_2nd) + return white_noise, count_cws_50, count_cws_01, hc_ss, hc_bg, dp_max, dp_2nd + +def get_noise_arrays_temp(target, red=False): + temp_name = '/Users/emigardiner/GWs/holodeck/output/temp' + temp_name += f'/noise_arrays_{target}.npz' + file = np.load(temp_name) + white_noise = file['white_noise'] + count_cws_50 = file['count_cws_50'] + count_cws_01 = file['count_cws_01'] + hc_ss = file['hc_ss'] + hc_bg = file['hc_bg'] + dp_max = file['dp_max'] + dp_2nd = file['dp_2nd'] + if red: + red_noise = file['red_noise'] + file.close() + return white_noise, red_noise, count_cws_50, count_cws_01, hc_ss, hc_bg, dp_max, dp_2nd + file.close() + return white_noise, count_cws_50, count_cws_01, hc_ss, hc_bg, dp_max, dp_2nd + + + +def build_anis_var_arrays( + target, nvars=21, nreals=500, shape=None, + gw_only=False, + nloudest=10, + lmax=8, nside=8, + var_hard_time=None, + figpath = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/anis_var', + +): + """ Calculate xx=params and yy=C1C0 + + """ + data, params, = get_data(target, nvars=nvars, nreals=nreals, shape=shape, + gw_only=gw_only, + nloudest=nloudest, dets=False, + var_hard_time=var_hard_time) + xx=[] + yy=[] + cl=[] + for pp, par in enumerate(params): + xx.append(params[pp][target]) + _, Cl = holo.anisotropy.sph_harm_from_hc( + data[pp]['hc_ss'], data[pp]['hc_bg'], nside=nside, lmax=lmax + ) + yy.append(Cl[...,1]/Cl[...,0]) + cl.append(Cl) + + filename = figpath+f'/anis_var_arrays_{target}' + filename += f"_l{lmax}_ns{nside}" + filename = append_filename( + filename, + gw_only=gw_only, + nloudest=nloudest, bgl=nloudest, + var_hard_time=var_hard_time) + filename += '.npz' + + np.savez(filename, xx_params=xx, yy_c1c0=yy, cl=cl) + return xx, yy, cl + +def get_anis_var_arrays( + target, + gw_only=False, + nloudest=10, + lmax=8, nside=8, + var_hard_time=None, + figpath = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/anis_var', + +): + """ Get xx=params and yy=C1C0 + + """ + + filename = figpath+f'/anis_var_arrays_{target}' + filename += f"_l{lmax}_ns{nside}" + filename = append_filename( + filename, + gw_only=gw_only, + nloudest=nloudest, bgl=nloudest, + var_hard_time=var_hard_time) + filename += '.npz' + + file = np.load(filename) + xx_params = file['xx_params'] + yy_c1c0 = file['yy_c1c0'] + cl = file['cl'] + file.close() + return xx_params, yy_c1c0, cl + + +def build_anis_freq_arrays( + target, nvars=21, nreals=500, shape=None, + gw_only=False, nloudest=10, + parvars = [0,10,20], + lmax=8, nside=8, + var_hard_time=None, + + ): + + if np.any(np.array(parvars)>nvars): + parvars = np.arange(nvars) + print(f'setting new parvars to {parvars}') + + data, params, = get_data(target, dets=False, + nvars=nvars, nreals=nreals, shape=shape, # keep as defaults + gw_only=gw_only, nloudest=nloudest, + var_hard_time=var_hard_time + ) + + + yy_cl = [] # PV len array of (F,R,l=1) arrays + params_cl = [] + xx_fobs = data[0]['fobs_cents'] + for var in parvars: + _, Cl = holo.anisotropy.sph_harm_from_hc( + data[var]['hc_ss'], data[var]['hc_bg'], nside=nside, lmax=lmax + ) + yy_cl.append(Cl) + params_cl.append(params[var]) + + filename = f'/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/anis_freq/anis_freq_{target}' + filename = append_filename(filename, nloudest=nloudest, bgl=nloudest, + gw_only=gw_only, + var_hard_time=var_hard_time) + filename += f"_pv{len(parvars)}" + if nside != 8: + filename += f"_ns{nside}" + if lmax != 8: + filename += f"_lmax{lmax}" + filename += f".npz" + np.savez(filename, yy_cl=yy_cl, xx_fobs=xx_fobs, params_cl=params_cl) + + return yy_cl, xx_fobs, params_cl + +def get_anis_freq_arrays( + target, nvars=21, nreals=500, nskies=100, shape=None, + gw_only=False, + nloudest=10, + parvars = [0,10,20], + nside=8, lmax=8, + var_hard_time=None, + + ): + + + if np.any(np.array(parvars)>nvars): + parvars = np.arange(nvars) + + filename = f'/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/anis_freq/anis_freq_{target}' + filename = append_filename(filename, nloudest=nloudest, bgl=nloudest, + gw_only=gw_only, + var_hard_time=var_hard_time) + filename += f"_pv{len(parvars)}" + if nside != 8: + filename += f"_ns{nside}" + if lmax != 8: + filename += f"_lmax{lmax}" + filename += f".npz" + file = np.load(filename, allow_pickle=True) + yy_cl=file['yy_cl'] + xx_fobs=file['xx_fobs'] + params_cl=file['params_cl'] + + return yy_cl, xx_fobs, params_cl + + +def build_hcpar_arrays( + target,nvars=21, nreals=500, shape=None, + gw_only=False, + nloudest=1, + parvars = [0,1,2,3,4], ss_zero=True, + var_hard_time=None, + ): + """ Save and return hcpar arrays for plotting + + returns + ----- + xx : + yy_ss : len(parvars) array of [F,R] NDarrays + 0th loudest [hc_ss, mass, distance] in dimensionless, M_sol, and Mpc units + yy_bg : len(parvars) array of [F,R] NDarrays + Background (all but nloudets) average [hc_ss, mass, distance] + in dimensionless, M_sol, and Mpc units + labels : nparvars array of target values + + """ + + if np.any(np.array(parvars)>nvars): + parvars = np.arange(nvars) + + labels = [] + yy_ss = [] + yy_bg = [] + data, params, = get_data(target, dets=False, + nvars=nvars, nreals=nreals, shape=shape, # keep as defaults + gw_only=gw_only, + nloudest=nloudest, + var_hard_time=var_hard_time) + + fobs_cents = data[0]['fobs_cents'] + xx = fobs_cents * YR + + for vv, var in enumerate(parvars): + labels.append(f"{params[var][target]}") + + hc_ss = data[var]['hc_ss'] + hc_bg = data[var]['hc_bg'] + + sspar = data[var]['sspar'] + bgpar = data[var]['bgpar'] + + sspar = holo.single_sources.all_sspars(fobs_cents, sspar) + bgpar = bgpar*holo.single_sources.par_units[:,np.newaxis,np.newaxis] + sspar = sspar*holo.single_sources.par_units[:,np.newaxis,np.newaxis,np.newaxis] + + # parameters to plot + if ss_zero: + _yy_ss = [hc_ss[...,0], sspar[0,...,0], #sspar[1,...,0], # sspar[2,], # strain, mass, mass ratio, + sspar[4,...,0]] # final comoving distance, single loudest only + else: + _yy_ss = [hc_ss[...], sspar[0,...], #sspar[1,...,0], # sspar[2,], # strain, mass, mass ratio, + sspar[4,...]] # final comoving distance, single loudest only + + _yy_bg = [hc_bg, bgpar[0], #bgpar[1], # strain, mass, mass ratio, initial redshift, final com distance + bgpar[4],] + yy_ss.append(_yy_ss) + yy_bg.append(_yy_bg) + + filename = f'/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/hcpar/hcpar_{target}' + filename = append_filename(filename, nloudest=nloudest, gw_only=gw_only, bgl=nloudest, + var_hard_time=var_hard_time) + filename += f"_pv{len(parvars)}" + filename += f".npz" + + np.savez(filename, xx=xx, yy_ss=yy_ss, yy_bg=yy_bg, labels=labels) + return np.array(xx), np.array(yy_ss), np.array(yy_bg), labels + +def get_hcpar_arrays( + target, + gw_only=False, + nloudest=1, + parvars = [0,1,2,3,4], + nvars=21, + var_hard_time=None, + ): + """ Save and return hcpar arrays for plotting + + returns + ------- + xx : [F,] 1Darray + yy_ss : len(parvars) array of [F,R] NDarrays + 0th loudest [hc_ss, mass, distance] in dimensionless, M_sol, and Mpc units + yy_bg : len(parvars) array of [F,R] NDarrays + Background (all but nloudets) average [hc_ss, mass, distance] + in dimensionless, M_sol, and Mpc units + labels : nparvars array of target values + + """ + + + if np.any(np.array(parvars)>nvars): + parvars = np.arange(nvars) + + filename = f'/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/hcpar/hcpar_{target}' + filename = append_filename(filename, nloudest=nloudest, gw_only=gw_only, bgl=nloudest, + var_hard_time=var_hard_time) + filename += f"_pv{len(parvars)}" + filename += f".npz" + + file = np.load(filename, allow_pickle=True) + xx = file['xx'] + yy_ss = file['yy_ss'] + yy_bg = file['yy_bg'] + labels = file['labels'] + + return xx, yy_ss, yy_bg, labels + +def build_favg_arrays( + target, nreals=500, nskies=100, + gw_only=False, red2white=None, red_gamma=None, + nloudest=10, bgl=10, cv=None, + gsc_flag=False, dsc_flag=False, divide_flag=False, nexcl=0, + var_hard_time=None, + figpath = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/favg',): + + data, params, dsdat = get_data(target, + gw_only=gw_only, red2white=red2white, red_gamma=red_gamma, + nloudest=nloudest, bgl=bgl, nreals=nreals, nskies=nskies, + gsc_flag=gsc_flag, dsc_flag=dsc_flag, divide_flag=divide_flag, nexcl=nexcl, + var_hard_time=var_hard_time) + + xx = [] # param + favg = [] # frequency means in log space + stdv = [] # stdev in log space + + freqs = data[0]['fobs_cents'] + nfreqs = len(freqs) + freqs = np.repeat(freqs, nreals*nskies*nloudest).reshape( + nfreqs, nreals, nskies, nloudest) + + for pp, par in enumerate(params): + xx.append(params[pp][target]) + dpssi = dsdat[pp]['gamma_ssi'] + logmean, logvar2 = weighted_mean_variance(np.log10(freqs), weights=dpssi) + + favg.append(logmean) + stdv.append(np.sqrt(logvar2)) + + filename = figpath+f'/favg_{target}' + filename = append_filename(filename, + gw_only=gw_only, red_gamma=red_gamma, red2white=red2white, + nloudest=nloudest, bgl=bgl, cv=cv, + gsc_flag=gsc_flag, dsc_flag=dsc_flag, divide_flag=divide_flag, nexcl=nexcl, + var_hard_time=var_hard_time) + + filename=filename+'.npz' + np.savez(filename, xx = xx, yy_log = favg, sd_log=stdv) + return xx, favg, stdv + +def get_favg_arrays( + target, + gw_only=False, red2white=None, red_gamma=None, + nloudest=10, bgl=10, cv=None, + gsc_flag=False, dsc_flag=False, divide_flag=False, nexcl=0, + var_hard_time=None, + figpath = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/favg',): + + + filename = figpath+f'/favg_{target}' + filename = append_filename(filename, + gw_only=gw_only, red_gamma=red_gamma, red2white=red2white, + nloudest=nloudest, bgl=bgl, cv=cv, + gsc_flag=gsc_flag, dsc_flag=dsc_flag, divide_flag=divide_flag, nexcl=nexcl, + var_hard_time=var_hard_time) + + filename=filename+'.npz' + file = np.load(filename) + xx = file['xx'] + yy_log = file['yy_log'] + sd_log = file['sd_log'] + file.close() + + return xx, yy_log, sd_log + + +def get_dp_arrays( + target, nvars=21, nreals=500, nskies=100, shape=None, debug=False, + red=False, cv='midclbrt', gw_only=False, bgl=1, + ): + + path = '/Users/emigardiner/GWs/holodeck/output/anatomy_redz/figdata/dpboth' + filename = path+f'/dp_arrays_{target}_cv{cv}' + if bgl != 10: + filename += f"_bgl{bgl}" + if gw_only: + filename = filename+'_gw' + filename = filename + '.npz' + + file = np.load(filename) + if debug: print(f"{filename}\n{file.files}") + xx = file['xx_params'] + yy_ss = file['yy_ss'] + ev_ss = file['ev_ss'] + yy_bg = file['yy_bg'] + file.close() + return xx, yy_ss, ev_ss, yy_bg \ No newline at end of file diff --git a/holodeck/extensions.py b/holodeck/extensions.py index 7c66a794..0081ac96 100644 --- a/holodeck/extensions.py +++ b/holodeck/extensions.py @@ -2,8 +2,11 @@ """ import holodeck as holo -from holodeck import log +from holodeck import log, cosmo, gravwaves from holodeck.constants import MSOL, GYR +import numpy as np +import kalepy as kale +from holodeck.sams import cyutils as sam_cyutils class Realizer: @@ -53,3 +56,226 @@ def __call__(self, down_sample=None): weights = self._binary_weights samples = evo._sample_universe__resample(fobs_orb_edges, vals, weights, down_sample) return names, samples + + +class Realizer_SAM: + def __init__( + self, fobs_orb_edges, sam=None, hard=None, params=None, + pspace=holo.param_spaces.PS_Uniform_09B(holo.log, nsamples=1, sam_shape=None, seed=None)): + """Construct a Realizer for a given semi-analytic model and hardening model, + or build this model using params and a pspace. + + Parameters + ---------- + sam : Semi_Analytic_Model object or None + Semi-analytic model instance, if not using pspace. + hard : Fixed_Time_2PL_SAM object, GW_Only object, None + Hardening model instance, if not using pspace. + params : dict or None + Parameters for a given parameter space, if sam is not provided. + pspace : _Param_Space object + Parameter space. + + + NOTE: To match the Realizer above I could initialize with weights and whatnot, then + possibly use the same resample/downsample function. + """ + + # check that ('sam' and 'hard') OR 'params' is provided + if params is not None: + if sam is not None or hard is not None: + err = "Only 'params' or ('sam' and 'hard') should be provided." + raise ValueError(err) + sam, hard = pspace.model_for_params(params=params, sam_shape=pspace.sam_shape,) + else: + if sam is None or hard is None: + err = "'params' or ('sam' and 'hard') must be provided." + raise ValueError(err) + + self._sam = sam + self._hard = hard + self._fobs_orb_edges = fobs_orb_edges + + def __call__(self, nreals=100, clean=False): + """ Calculate samples and weights for an entire semi-analytic population. + + Parameters + ---------- + nreals : int + Number of realizations + clean : boolean + Whether or not to make a samples array for every realization + and clean weights==zero bins from each array + + Returns + ------- + names : array of strings + Names of the parameters returned in samples + samples : array of R or 4 NDarrays + if clean: R arrays of 4 x N_clean NDarrays [R,] x [4,N_clean] for each realization + else: NDarrays for mass, ratio, redshift, and frequency [4,M*Q*Z*F] + weights : array of R NDarrays + array of number of sources per sample bin for R + If clean, the shape is [R,] arrays of len N_clean for each realizations, with zero values removed. + Otherwise, the shape is [R, M*Q*Z*F] + + """ + + sam = self._sam + hard = self._hard + fobs_orb_edges = self._fobs_orb_edges + + fobs_orb_cents = kale.utils.midpoints(fobs_orb_edges) + fobs = 2.0 * fobs_orb_cents + + + # ---- Calculate number of binaries in each bin + + redz, diff_num = sam_cyutils.dynamic_binary_number_at_fobs( + fobs_orb_cents, sam, hard, cosmo + ) + + edges = [sam.mtot, sam.mrat, sam.redz, fobs_orb_edges] + number = sam_cyutils.integrate_differential_number_3dx1d(edges, diff_num) # fractional number per bin + + samples = get_samples_from_edges(edges, redz, number.shape, flatten=True) + names = ['mtot', 'mrat', 'redz', 'fobs'] + number = number.flatten() + shape = (number.size, nreals) + weights = gravwaves.poisson_as_needed(number[..., np.newaxis] * np.ones(shape)).reshape(shape) + + + + if clean: + nonzero_samples = [] + nonzero_weights = [] + for rr in range(nreals): + nonzero = weights[:,rr]!=0 + mtot = samples[0][nonzero] + mrat = samples[1][nonzero] + redz = samples[2][nonzero] + fobs = samples[3][nonzero] + nonzero_samples.append([mtot, mrat, redz, fobs]) + nonzero_weights.append(weights[:,rr][nonzero]) + + weights = nonzero_weights + samples = nonzero_samples + + return names, samples, weights + + +def get_samples_from_edges(edges, redz, number_shape, flatten=True): + """ Get the sample parameters for every bin center and return as flattened arrays. + + Parameters + ---------- + edges : array of [M+1,], [Q+1,], [Z+1,], and [F+1,] NDarrays + Edges for mtot, mrat, redz, fobs_orb_edges + redz : [M+1, Q+1, Z+1, F+1] NDarray + Final redshifts + number_shape : array + Shape [M,Q,Z,F] + flatten : boolean + Whether or not to flatten each sample array + + Returns + ------- + samples : array of 4 flattened [M*Q*Z*F,] NDarrays + + """ + + # ---- Find bin center properties + mtot = kale.utils.midpoints(edges[0]) #: total mass + mrat = kale.utils.midpoints(edges[1]) #: mass ratio + fobs_orb_cents = kale.utils.midpoints(edges[3]) + fobs = 2.0 * fobs_orb_cents #: gw fobs + + for dd in range(3): + redz = np.moveaxis(redz, dd, 0) + redz = kale.utils.midpoints(redz, axis=0) # get final redz at bin centers + redz = np.moveaxis(redz, 0, dd) + sel = (redz > 0.0) # identify emitting sources + redz[~sel] = -1.0 # set all other redshifts to zero + redz[redz<0] = -1.0 + + # get bin shape + nmtot = len(mtot) + nmrat = len(mrat) + nredz = redz.shape[2] + nfobs = len(fobs) + # if np.any([number_shape[0] != nmtot, number_shape[1]!=nmrat, number_shape[2]!=nredz, number_shape[3]!=nfobs]): + # err = f"Parameter bin shape [{nmtot=}, {nmrat=}, {nredz=}, {nfobs=}] does not match {number_shape=}." + # raise ValueError(err) + print(f"Parameter bin shape [{nmtot=}, {nmrat=}, {nredz=}, {nfobs=}] should match {number_shape=}.") + + # Reshape arrays to [M,Q,Z,F] + mtot = np.repeat(mtot, nmrat*nredz*nfobs).reshape(nmtot, nmrat, nredz, nfobs) + + mrat = np.repeat(mrat, nmtot*nredz*nfobs).reshape(nmrat, nmtot, nredz, nfobs) # Q,M,Z,F + mrat = np.swapaxes(mrat, 0, 1) # M,Q,Z,F + + fobs = np.repeat(fobs, nmrat*nredz*nmtot).reshape(nfobs, nmrat, nredz, nmtot) # F,Q,Z,M + fobs = np.swapaxes(fobs, 0, 3) # M,Q,Z,F + + # check shapes again + if np.any([mtot.shape != number_shape, + mrat.shape != number_shape, + redz.shape != number_shape, + fobs.shape != number_shape]): + err = f"Sample shapes don't all match number! {mtot.shape=}, {mrat.shape=}, {redz.shape=}, {fobs.shape=}" + raise ValueError(err) + + if flatten: + samples = [mtot.flatten(), mrat.flatten(), redz.flatten(), fobs.flatten()] + else: + samples = [mtot, mrat, redz, fobs] + + return samples + + +def realizer_single_sources(params, nreals, nloudest, nfreqs=40, log10=False, + pspace = holo.param_spaces.PS_Uniform_09B(holo.log, nsamples=1, sam_shape=None, seed=None)): + """ Like Realizer but using single sources from a SAM instead of Illustris populations + + Parameters + ---------- + params : dict + model parameters for parameter space + nreal : int + number of realizations + nloudest : int + number of loudest sources to use + + Returns + ------- + names : names of parameters + samples : [R, 4, nloudest*nreals] NDarray + mtot, mrat, redz, and fobs of each source in log space + + sspar is in shape [F,R,L] + + """ + fobs_cents, fobs_edges = holo.utils.pta_freqs(num=nfreqs) + + sam, hard = pspace.model_for_params(params=params, sam_shape=None,) + _, _, sspar, bgpar = sam.gwb( + fobs_edges, hard=hard, nreals=nreals, nloudest=nloudest, params=True) + + vals_names=['mtot', 'mrat', 'redz', 'fobs'] + fobs = np.repeat(fobs_cents, nreals*nloudest).reshape(nfreqs, nreals, nloudest) + mtot = sspar[0] # g + mrat = sspar[1] + redz = sspar[3] # final redshift, not initial + + samples = [] + for par in [mtot, mrat, redz, fobs]: # starts in shape F,R,L + par = np.swapaxes(par, 0, 1) # R,F,L + par = par.reshape(nreals, nfreqs*nloudest) # R, F*L + if log10: + samples.append(np.log10(par)) + else: + samples.append(par) + samples = np.array(samples) # 4,R,F*L + samples = np.swapaxes(samples, 0,1) # R,4,F*L + + return vals_names, samples \ No newline at end of file diff --git a/holodeck/plot.py b/holodeck/plot.py index 252db98a..66474722 100644 --- a/holodeck/plot.py +++ b/holodeck/plot.py @@ -38,7 +38,7 @@ LABEL_SEPARATION_PC = r"Binary Separation $[\mathrm{pc}]$" LABEL_CHARACTERISTIC_STRAIN = r"GW Characteristic Strain" LABEL_HARDENING_TIME = r"Hardening Time $[\mathrm{Gyr}]$" - +LABEL_CLC0 = r"$C_\ell / C_0$" PARAM_KEYS = { 'hard_time': r"phenom $\tau_f$", @@ -58,7 +58,8 @@ } LABEL_DPRATIO = r"$\langle N_\mathrm{SS} \rangle / \mathrm{DP}_\mathrm{BG}$" -LABEL_EVSS = r"$\langle N_\mathrm{SS} \rangle" +LABEL_EVSS = r"$\langle N_\mathrm{SS} \rangle$" +LABEL_DPBG = r"$\mathrm{DP}_\mathrm{BG}$" COLORS_MPL = plt.rcParams['axes.prop_cycle'].by_key()['color'] @@ -752,7 +753,8 @@ def draw_med_conf(ax, xx, vals, fracs=[0.50, 0.90], weights=None, plot={}, fill= return (hh, gg) -def draw_med_conf_color(ax, xx, vals, fracs=[0.50, 0.90], weights=None, plot={}, fill={}, filter=False, color=None): +def draw_med_conf_color(ax, xx, vals, fracs=[0.50, 0.90], weights=None, plot={}, fill={}, + filter=False, color=None, linestyle='-'): plot.setdefault('alpha', 0.75) fill.setdefault('alpha', 0.2) percs = np.atleast_1d(fracs) @@ -776,7 +778,7 @@ def draw_med_conf_color(ax, xx, vals, fracs=[0.50, 0.90], weights=None, plot={}, # plot median if color is not None: - hh, = ax.plot(xx, med, color=color, **plot) + hh, = ax.plot(xx, med, color=color, linestyle=linestyle, **plot) else: hh, = ax.plot(xx, med, **plot) @@ -1459,6 +1461,17 @@ def _contour2d(ax, edges, hist, levels, outline=True, **kwargs): return edges, hist, cont +def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100): + ''' + https://stackoverflow.com/a/18926541 + ''' + if isinstance(cmap, str): + cmap = plt.get_cmap(cmap) + new_cmap = mpl.colors.LinearSegmentedColormap.from_list( + 'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval), + cmap(np.linspace(minval, maxval, n))) + return new_cmap + # ================================================================================================= # ==== Below Needs Review / Cleaning ==== # ================================================================================================= diff --git a/holodeck/sams/cyutils.pyx b/holodeck/sams/cyutils.pyx index 28976d3f..1ef851db 100644 --- a/holodeck/sams/cyutils.pyx +++ b/holodeck/sams/cyutils.pyx @@ -762,195 +762,3 @@ cdef int _dynamic_binary_number_at_fobs_gw( return 0 - - -# ================================================================================================== -# ==== DetStats Functions ==== -# ================================================================================================== - - -def gamma_of_rho_interp(rho, rsort, rho_interp_grid, gamma_interp_grid): - """ - rho : 1Darray of scalars - SNR of single sources, in flat array - rsort : 1Darray - order of flat rho values smallest to largest - rho_interp_grid : 1Darray - rho values corresponding to each gamma - gamma_interp_grid : 1Darray - gamma values corresponding to each rho - - """ - # pass in the interp grid - cdef np.ndarray[np.double_t, ndim=1] gamma = np.zeros(rho.shape) - - _gamma_of_rho_interp(rho, rsort, rho_interp_grid, gamma_interp_grid, gamma) - - return gamma - -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.nonecheck(False) -@cython.cdivision(True) -cdef int _gamma_of_rho_interp( - double[:] rho, long[:] rsort, - double[:] rho_interp_grid, double[:] gamma_interp_grid, - # output - double[:] gamma - ): - """ Find gamma of rho by interpolation over rho and gamma grids. - """ - - cdef int n_rho = rho.size - cdef int n_interp = rho_interp_grid.size - cdef int ii, kk, rr - ii = 0 # get rho in order using rho[rsort[ii]] - - for kk in range(n_rho): - rr = rsort[kk] # index of next largest rho, equiv to rev in redz calculation - # print('kk =',kk,' rr =', rr, 'rho[rr] =', rho[rr]) - # get to the right index of the interpolation-grid - while (rho_interp_grid[ii+1] < rho[rr]) and (ii < n_interp -2): - ii += 1 - # print('ii =',ii, ' rho_interp[ii] =', rho_interp_grid[ii], ' rho_interp[ii+1] =', rho_interp_grid[ii+1]) - # interpolate - gamma[rr] = interp_at_index(ii, rho[rr], rho_interp_grid, gamma_interp_grid) - # print('rho =', rho[rr], ' gamma =', gamma[rr], '\n')\ - - return 0 - - -def snr_ss(amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs): - """ Calculate single source SNR - - - Parameters - ---------- - amp : (F,R,L) NDarray - Dimensionless strain amplitude of loudest single sources - F_iplus : (P,F,S,L) NDarray - Antenna pattern function for each pulsar. - F_icross : (P,F,S,L) NDarray - Antenna pattern function for each pulsar. - iotas : (F,S,L) NDarray - Inclination, used to calculate: - a_pol = 1 + np.cos(iotas) **2 - b_pol = -2 * np.cos(iotas) - dur : scalar - Duration of observations. - Phi_0 : (F,S,L) NDarray - Initial GW phase - S_i : (P,F,R,L) NDarray - Total noise of each pulsar wrt detection of each single source, in s^3 - freqs : (F,) 1Darray - Observed frequency bin centers. - - Returns - ------- - snr_ss : (F,R,S,L) NDarray - SNR from the whole PTA for each single source with - each realized sky position (S) and realized strain (R) - - """ - nfreqs, nreals, nloudest = amp.shape[0], amp.shape[1], amp.shape[2] - npsrs, nskies = F_iplus.shape[0], F_iplus.shape[2] - cdef np.ndarray[np.double_t, ndim=4] snr_ss = np.zeros((nfreqs, nreals, nskies, nloudest)) - _snr_ss( - amp, F_iplus, F_icross, iotas, dur, Phi_0, S_i, freqs, - npsrs, nfreqs, nreals, nskies, nloudest, - snr_ss) - return snr_ss - -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.nonecheck(False) -@cython.cdivision(True) -cdef int _snr_ss( - double[:,:,:] amp, - double[:,:,:,:] F_iplus, - double[:,:,:,:] F_icross, - double[:,:,:] iotas, - double dur, - double[:,:,:] Phi_0, - double[:,:,:,:] S_i, - double[:] freqs, - long npsrs, long nfreqs, long nreals, long nskies, long nloudest, - # output - double[:,:,:,:] snr_ss - ): - """ - - Parameters - ---------- - amp : (F,R,L) NDarray - Dimensionless strain amplitude of loudest single sources - F_iplus : (P,F,S,L) NDarray - Antenna pattern function for each pulsar. - F_icross : (P,F,S,L) NDarray - Antenna pattern function for each pulsar. - iotas : (F,S,L) NDarray - Inclination, used to calculate: - a_pol = 1 + np.cos(iotas) **2 - b_pol = -2 * np.cos(iotas) - dur : scalar - Duration of observations. - Phi_0 : (F,S,L) NDarray - Initial GW phase - S_i : (P,F,R,L) NDarray - Total noise of each pulsar wrt detection of each single source, in s^3 - freqs : (F,) 1Darray - Observed frequency bin centers. - snr_ss : (F,R,S,L) NDarray - Pointer to single source SNR array, to be calculated. - - NOTE: This may be improved by moving some of the math outside the function. - I.e., passing in sin/cos of NDarrays to be used. - """ - - cdef int pp, ff, rr, ss, ll - cdef float a_pol, b_pol, Phi_T, pta_snr_sq, coef, term1, term2, term3 - # print('npsrs %d, nfreqs %d, nreals %d, nskies %d, nloudest %d' % (npsrs, nfreqs, nreals, nskies, nloudest)) - - for ff in range(nfreqs): - for ss in range(nskies): - for ll in range(nloudest): - a_pol = 1 + pow(cos(iotas[ff,ss,ll]), 2.0) - b_pol = -2 * cos(iotas[ff,ss,ll]) - Phi_T = 2 * M_PI * freqs[ff] * dur + Phi_0[ff,ss,ll] - for rr in range(nreals): - pta_snr_sq = 0 - for pp in range(npsrs): - # calculate coefficient depending on - # function of amp, S_i, and freqs - coef = pow(amp[ff,rr,ll], 2.0) / (S_i[pp,ff,rr,ll] * 8 * pow(M_PI * freqs[ff], 3.0)) - - # calculate terms that depend on p, f, s, and l - # functions of F_iplus, F_icross, a_pol, b_pol, Phi_0, and Phi_T - term1 = ( - pow(a_pol * F_iplus[pp,ff,ss,ll], 2.0) - * (Phi_T * (1.0 + 2.0 * pow(sin(Phi_0[ff,ss,ll]), 2.0)) - + cos(Phi_T) * (-1.0 * sin(Phi_T) + 4.0 * cos(Phi_0[ff,ss,ll])) - - 4.0 * sin(Phi_0[ff,ss,ll]) - ) - ) - term2 = ( - pow(b_pol * F_icross[pp,ff,ss,ll], 2.0) - * (Phi_T * (1.0 + 2.0 * pow(cos(Phi_0[ff,ss,ll]), 2.0)) - + sin(Phi_T) * cos(Phi_T) - 4.0 * cos(Phi_0[ff,ss,ll]) - ) - ) - term3 = ( - -2.0 * a_pol * b_pol * F_iplus[pp,ff,ss,ll] * F_icross[pp,ff,ss,ll] - * (2.0 * Phi_T * sin(Phi_T) *cos(Phi_0[ff,ss,ll]) - + sin(Phi_T) * (sin(Phi_T) - 2.0 * sin(Phi_0[ff,ss,ll]) - + 2.0 * cos(Phi_T) * cos(Phi_0[ff,ss,ll]) - - 2.0 * cos(Phi_0[ff,ss,ll]) - ) - ) - ) - pta_snr_sq += coef*(term1 + term2 + term3) # sum snr^2 of all pulsars for a single source - - # set snr for a single source, using sum from all pulsars - snr_ss[ff,rr,ss,ll] = sqrt(pta_snr_sq) - - diff --git a/holodeck/single_sources.py b/holodeck/single_sources.py index a0287c19..07195290 100644 --- a/holodeck/single_sources.py +++ b/holodeck/single_sources.py @@ -96,8 +96,6 @@ def ss_gws_redz(edges, redz, number, realize, loudest = 1, params = False): err = np.sum(np.logical_and(redz<0, redz!=-1)) err = f"{err} redz < 0 and !=-1 found in redz, in ss_gws_redz()" raise ValueError(err) - - print('passed redz check at the beginning of ss_gws_redz()') # For multiple realizations, using cython if(utils.isinteger(realize)): @@ -2141,46 +2139,21 @@ def plot_params(axs, xx, REALS=1, LABEL='', grid=None, if(SHOW_LEGEND): axs[ii,jj].legend(loc='lower left') - ################################################### -############ DETECTION STATISTICS ################# +############ Utilities ################# ################################################### +def resample_loudest(hc_ss, hc_bg, nloudest): + if nloudest > hc_ss.shape[-1]: # check for valid nloudest + err = f"{nloudest=} for detstats must be <= nloudest of hc data" + raise ValueError(err) + + # recalculate new hc_bg and hc_ss + new_hc_bg = np.sqrt(hc_bg**2 + np.sum(hc_ss[...,nloudest:-1]**2, axis=-1)) + new_hc_ss = hc_ss[...,0:nloudest] -def threshold_hc(): - """ - Rosado+ 2015, SNR calculation - S := cross correlation between pulsars - S = \int_{-T/2}^{T/2} dt \int dt' s_i(t) s_j(t') Q(t,t') - where T is the observation time, s_i(t) and s_j(t) are the data - from two different pulsars, Q(t,t') is a filter function. + return new_hc_ss, new_hc_bg - S_T := pre-defined detection threshold - """ - return 0 -def ss_occurence_rate(hc_ss, S_T): - """ - Parameters: - ------------ - hc_ss : (F, R) NDarray of scalars - characteristic strain of single sources - S_T : float - threshold strain? - - """ - return 0 - -def bg_occurence_rate(): - """ - According to Rosado+ 2015 - Universe may contain GWB if S >= S_T - """ - -def false_alarm_probability(): - """ - TODO - """ - return 0 \ No newline at end of file