diff --git a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py index f5c28378..b7a92aff 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/lalsimutils.py @@ -33,6 +33,7 @@ import EOBRun_module has_external_teobresum=True except: + #print('Failed to import EOBRun_module') has_external_teobresum=False True; # print(" - no EOBRun (TEOBResumS) - ") info_use_resum_polarizations = False @@ -331,7 +332,7 @@ def lsu_StringFromPNOrder(order): # # Class to hold arguments of ChooseWaveform functions # -valid_params = ['m1', 'm2', 's1x', 's1y', 's1z', 's2x', 's2y', 's2z', 'chi1_perp', 'chi2_perp', 'chi1_perp_bar', 'chi2_perp_bar','chi1_perp_u', 'chi2_perp_u', 's1z_bar', 's2z_bar', 'lambda1', 'lambda2', 'theta','phi', 'phiref', 'psi', 'incl', 'tref', 'dist', 'mc', 'mc_ecc', 'eta', 'delta_mc', 'chi1', 'chi2', 'thetaJN', 'phiJL', 'theta1', 'theta2', 'cos_theta1', 'cos_theta2', 'theta1_Jfix', 'theta2_Jfix', 'psiJ', 'beta', 'cos_beta', 'sin_phiJL', 'cos_phiJL', 'phi12', 'phi1', 'phi2', 'LambdaTilde', 'DeltaLambdaTilde', 'lambda_plus', 'lambda_minus', 'q', 'mtot','xi','chiz_plus', 'chiz_minus', 'chieff_aligned','fmin','fref', "SOverM2_perp", "SOverM2_L", "DeltaOverM2_perp", "DeltaOverM2_L", "shu","ampO", "phaseO",'eccentricity','chi_pavg','mu1','mu2','eos_table_index', 'E0', 'p_phi0'] +valid_params = ['m1', 'm2', 's1x', 's1y', 's1z', 's2x', 's2y', 's2z', 'chi1_perp', 'chi2_perp', 'chi1_perp_bar', 'chi2_perp_bar','chi1_perp_u', 'chi2_perp_u', 's1z_bar', 's2z_bar', 'lambda1', 'lambda2', 'theta','phi', 'phiref', 'psi', 'incl', 'tref', 'dist', 'mc', 'mc_ecc', 'eta', 'delta_mc', 'chi1', 'chi2', 'thetaJN', 'phiJL', 'theta1', 'theta2', 'cos_theta1', 'cos_theta2', 'theta1_Jfix', 'theta2_Jfix', 'psiJ', 'beta', 'cos_beta', 'sin_phiJL', 'cos_phiJL', 'phi12', 'phi1', 'phi2', 'LambdaTilde', 'DeltaLambdaTilde', 'lambda_plus', 'lambda_minus', 'q', 'mtot','xi','chiz_plus', 'chiz_minus', 'chieff_aligned','fmin','fref', "SOverM2_perp", "SOverM2_L", "DeltaOverM2_perp", "DeltaOverM2_L", "shu","ampO", "phaseO",'eccentricity','chi_pavg','mu1','mu2','eos_table_index', 'E0', 'p_phi0', 'hypclass', 'b_hyp', 'phi_scattter'] tex_dictionary = { "mtot": '$M$', @@ -374,6 +375,8 @@ def lsu_StringFromPNOrder(order): "eccentricity":"$e$", "E0": "$E_0 / M$", "p_phi0": r"$p_{\phi}^0$", + "b_hyp" : r"$b_{hyp}$", + "phi_scatter" : r"$\phi_{sc}$", # tex labels for inherited LI names "a1z": r'$\chi_{1,z}$', "a2z": r'$\chi_{2,z}$', @@ -884,7 +887,195 @@ def extract_param(self,p): if p == 'mtot': return (self.m2+self.m1) if p == 'q': - return self.m2/self.m1 + return self.m2/self.m1 + if p == 'hypclass': + # Checks type of hyperbolic waveform: scatter, plunge, zoomwhirl, or meaningless + # check if valid + if self.E0 == 0.0: + print('Invalid use of hypclass: non-hyperbolic configuration') + return None + + # Generate waveform + pars = { + 'M' : (self.m1+self.m2)/lal.MSUN_SI, + 'q' : self.m1/self.m2, + 'H_hyp' : self.E0, # energy at initial separation + 'j_hyp' : self.p_phi0, # angular momentum at initial separation + 'r_hyp' : 6000.0, + 'LambdaAl2' : self.lambda1, + 'LambdaBl2' : self.lambda2, + 'chi1' : self.s1z, #note that there are no transverse spins + 'chi2' : self.s2z, + 'dt' : self.deltaT, + 'domain' : 0, # 0 sets time domain + 'arg_out' : 1, # Request multipoles and dynamics as output of the function call - 1=yes + 'nqc' : 2, # sets the NQCs, 2=no + 'nqc_coefs_hlm' : 0, # Option for the NQC model used in the waveform. 0=none + 'nqc_coefs_flx' : 0, # Option for the NQC model used in the flux. 0=none + 'use_mode_lm' : [1], # 22 mode + 'output_lm' : [1], + 'srate_interp' : 1./self.deltaT, + 'use_geometric_units': 0, + 'interp_uniform_grid': 1, + 'initial_frequency' : self.fmin, + 'ode_tmax' : 3e4, + 'distance' : self.dist/(lal.PC_SI*1e6), + 'inclination' : self.incl, + 'output_hpc' : 0 # output plus and cross polarizations, 0=no + } + + t, hptmp, hctmp, hlmtmp, dym = EOBRun_module.EOBRunPy(pars) + + # peak finding to classifiy + + # wf amplitude - 22 mode only + tmp_22 = np.array(hlmtmp['1']) + distance_s = self.dist/lal.C_SI + m_total_s = MsunInSec*(self.m1+self.m2)/lal.MSUN_SI + M1=self.m1/lal.MSUN_SI + M2=self.m2/lal.MSUN_SI + nu=M1*M2/((M1+M2)**2) + tmp_22[0] *= (m_total_s/distance_s)*nu + amp = np.abs(tmp_22[0] * np.exp(-1j*(2*(np.pi/2.)+tmp_22[1]))) + amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding + + # Check for cases where the amplitude is max at the start or end + if np.argmax(amp_norm) == 0 or np.argmax(amp_norm) == len(amp_norm) - 1: + reclassify = True + else: + reclassify = False + + # peak finding to determine system type + height_thresh = 0.25 + prom_thresh = 0.1 + peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh) + peak_heights = props['peak_heights'] + # filtering out peaks so we only keep the local maxima + indices_to_keep = set() + sorted_indices = np.argsort(peak_heights)[::-1] + tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks. + for i in sorted_indices: + peak = peaks[i] + keep = True + for kept_index in indices_to_keep: + if abs(peaks[kept_index] - peak) <= tol: + keep = False + break + if keep: + indices_to_keep.add(i) + filtered_peaks = peaks[list(indices_to_keep)] + + # parsing number of peaks after filtering against distance tolerance + if len(filtered_peaks) == 1: + if np.abs(amp)[-1] > 1e-26: + # scatter waveform + return 'scatter' + else: + # plunge waveform + return 'plunge' + + elif len(filtered_peaks) == 0: + # meaningless waveform + reclassify = True + + else: + # zoom whirl waveform + return 'zoomwhirl' + + if reclassify: + print('Running re-classification') + # run minimal threshold peak finder + all_peaks, all_props = signal.find_peaks(amp_norm, height=0.0001, prominence=0.0001) + + # checking for minima if no peaks found + if len(all_peaks) == 0: + print("No peaks found, checking for minima instead...") + all_peaks, all_props = signal.find_peaks((-1*amp_norm + 1.0), height=0.0001, prominence=0.0001) + + if len(all_props['prominences']) > 3: + print('MANY peaks detected on reclassification, evaluating...') + # these can be scatter or plunge + if np.abs(amp)[-1] < 1e-26: + print('Reclassifying to Plunge') + return 'plunge' + else: + print('Reclassifying to Scatter') + return 'scatter' + elif len(all_props['prominences']) == 3 or len(all_props['prominences']) == 2 or len(all_props['prominences']) == 1: + # these are always scaters + print('Reclassifying to Scatter') + return 'scatter' + else: + print('No peaks detected after reclassifcation') + return 'meaningless' + + if p == 'b_hyp': + # Calculates impact parameter for hyperbolic cases. + # See https://arxiv.org/pdf/2309.07228 for details. + # check if valid use + if self.E0 == 0.0: + print('Invalid use of b_hyp: non-hyperbolic configuration') + return None + + h = self.E0 # mass normalized local system energy at initial separation + j = self.p_phi0 # mass normalized local system angular momentum at initial separation + + M1=self.m1/lal.MSUN_SI + M2=self.m2/lal.MSUN_SI + nu=M1*M2/((M1+M2)**2) + + E_eff = 1. + (h**2 - 1)/(2*nu) # effective energy of the one-body particle + b = j*h / (np.sqrt(E_eff**2 - 1)) # dimensionless impact parameter normalized by total mass + + return b + + if p == 'phi_scatter': + # Calculates scattering angle for hyperbolic scatters + # See https://arxiv.org/pdf/1402.7307, + # https://teobresums.bitbucket.io/gallery/scattering-angles/ + + # NOTE that you can calculate this for capture/plunge cases, but it is physically meaningless + + # check if valid use + if self.E0 == 0.0: + print('Invalid use of phi_scatter: non-hyperbolic configuration') + return None + + # Generate waveform + pars = { + 'M' : (self.m1+self.m2)/lal.MSUN_SI, + 'q' : self.m1/self.m2, + 'H_hyp' : self.E0, # energy at initial separation + 'j_hyp' : self.p_phi0, # angular momentum at initial separation + 'r_hyp' : 6000.0, + 'LambdaAl2' : self.lambda1, + 'LambdaBl2' : self.lambda2, + 'chi1' : self.s1z, #note that there are no transverse spins + 'chi2' : self.s2z, + 'dt' : self.deltaT, + 'domain' : 0, # 0 sets time domain + 'arg_out' : 1, # Request multipoles and dynamics as output of the function call - 1=yes + 'nqc' : 2, # sets the NQCs, 2=no + 'nqc_coefs_hlm' : 0, # Option for the NQC model used in the waveform. 0=none + 'nqc_coefs_flx' : 0, # Option for the NQC model used in the flux. 0=none + 'use_mode_lm' : [1], # 22 mode + 'output_lm' : [1], + 'srate_interp' : 1./self.deltaT, + 'use_geometric_units': 0, + 'interp_uniform_grid': 1, + 'initial_frequency' : self.fmin, + 'ode_tmax' : 3e4, + 'distance' : self.dist/(lal.PC_SI*1e6), + 'inclination' : self.incl, + 'output_hpc' : 0 # output plus and cross polarizations, 0=no + } + + t, hptmp, hctmp, hlmtmp, dyn = EOBRun_module.EOBRunPy(pars) + + angle = dyn['phi'][-1] - dyn['phi'][0] - np.pi + + return angle + if p == 'delta' or p=='delta_mc': # Same access routine return (self.m1-self.m2)/(self.m1+self.m2) if p == 'mc': @@ -1852,14 +2043,16 @@ def copy_lsctables_sim_inspiral(self, row): # Call the function to read lalmetaio.SimInspiral format self.copy_sim_inspiral(swigrow) - def scale_to_snr(self,new_SNR,psd, ifo_list,analyticPSD_Q=True): + def scale_to_snr(self,new_SNR,psd, ifo_list,analyticPSD_Q=True, **kwargs): """ scale_to_snr - evaluates network SNR in the ifo list provided (assuming *constant* psd for all..may change) - uses network SNR to rescale the distance of the source, so the SNR is now new_SNR - returns current_SNR, for sanity """ - deltaF=findDeltaF(self) + Lmax = kwargs.get('Lmax', 4) # Default to 4 if not specified in kwargs + + deltaF=findDeltaF(self, Lmax=Lmax) det_orig = self.detector IP = Overlap(fLow=self.fmin, fNyq=1./self.deltaT/2., deltaF=deltaF, psd=psd, full_output=True,analyticPSD_Q=analyticPSD_Q) @@ -1868,7 +2061,7 @@ def scale_to_snr(self,new_SNR,psd, ifo_list,analyticPSD_Q=True): for det in ifo_list: self.detector = det self.radec = True - h=hoff(self) + h=hoff(self, Lmax=Lmax) rho_ifo[det] = IP.norm(h) current_SNR_squared +=rho_ifo[det]*rho_ifo[det] current_SNR = np.sqrt(current_SNR_squared) @@ -2749,7 +2942,7 @@ def nextPow2(length): """ return int(2**np.ceil(np.log2(length))) -def findDeltaF(P): +def findDeltaF(P, **kwargs): """ Given ChooseWaveformParams P, generate the TD waveform, round the length to the next power of 2, @@ -2757,7 +2950,9 @@ def findDeltaF(P): This is useful b/c deltaF is needed to define an inner product which is needed for norm_hoft and norm_hoff functions """ - h = hoft(P) + Lmax = kwargs.get('Lmax', 4) # Default to 4 if not specified in kwargs + + h = hoft(P, Lmax=Lmax) return 1./(nextPow2(h.data.length) * P.deltaT) def estimateWaveformDuration(P,LmaxEff=2): @@ -2839,7 +3034,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs): extra_waveform_args.update(kwargs['extra_waveform_args']) extra_params = P.to_lal_dict_extended(extra_args_dict=extra_waveform_args) if P.approx==lalsim.TEOBResumS and has_external_teobresum and info_use_ext: - Lmax=4 + Lmax = kwargs.get('Lmax', 4) # Default to 4 if not specified in kwargs modes_used = [] distance_s = P.dist/lal.C_SI m_total_s = MsunInSec*(P.m1+P.m2)/lal.MSUN_SI @@ -2851,6 +3046,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs): M1=P.m1/lal.MSUN_SI M2=P.m2/lal.MSUN_SI nu=M1*M2/((M1+M2)**2) + hyp_wav = False print(P.eccentricity, P.E0) if (P.eccentricity == 0.0 and P.E0 == 0.0): print("Using ResumS master; not eccentric") @@ -2878,6 +3074,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs): } elif (P.eccentricity == 0.0): print("Using hyperbolic call RIFT O4b branch") + hyp_wav = True # convenient way to know if the waveform is hyperbolic pars = { 'M' : M1+M2, 'q' : M1/M2, @@ -2932,7 +3129,47 @@ def hoft(P, Fp=None, Fc=None,**kwargs): print("Starting EOBRun_module") t, hptmp, hctmp, hlmtmp, dyn = EOBRun_module.EOBRunPy(pars) print("EOBRun_module done") - hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2) + + if not(hyp_wav): + ## Set the epoch for non-hyperbolic cases ## + hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2) + else: + ## custom epoch for the hyperbolic case ## + # wf amplitude + amp = np.sqrt(np.abs(hptmp)**2+np.abs(hctmp)**2) + amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding + # peak finding to determine system type + height_thresh = 0.25*np.abs(amp_norm) + prom_thresh = 0.1*np.abs(amp_norm) + peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh) + peak_heights = props['peak_heights'] + # filtering out peaks so we only keep the local maxima + indices_to_keep = set() + sorted_indices = np.argsort(peak_heights)[::-1] + tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks. + for i in sorted_indices: + peak = peaks[i] + keep = True + for kept_index in indices_to_keep: + if abs(peaks[kept_index] - peak) <= tol: + keep = False + break + if keep: + indices_to_keep.add(i) + filtered_peaks = peaks[list(indices_to_keep)] + # parsing number of peaks after filtering against distance tolerance + if len(filtered_peaks) == 1: + # scatter case OR plunge case, we can set the epoch normally + hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2) + elif len(filtered_peaks) == 0: + # meaningless waveform, essentially a non-interacting case. Set the epoch normally + # These points should return very low likelihood and not interfere with the analysis + print('WARNING: no peak detected; non-interacting hyperbolic case') + hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2) + else: + # capture case, we need to force the epoch to be the last peak + hpepoch = -P.deltaT*filtered_peaks[-1] + hplen = len(hptmp) hp = {} hc = {} @@ -3013,15 +3250,56 @@ def hoft(P, Fp=None, Fc=None,**kwargs): if np.abs(value-ht.data.data[0]) > 0.01 * np.abs(ht.data.data[0]): n_samp=int(count/2) break - vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp))) + + # peak finding to determine system type + ht_norm = ht.data.data/np.amax(ht.data.data) + height_thresh = 0.25 + prom_thresh = 0.1 + + peaks, props = signal.find_peaks(ht_norm, height = height_thresh, prominence = prom_thresh) + + peak_heights = props['peak_heights'] + indices_to_keep = set() + sorted_indices = np.argsort(peak_heights)[::-1] + tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks. + for i in sorted_indices: + peak = peaks[i] + keep = True + for kept_index in indices_to_keep: + if abs(peaks[kept_index] - peak) <= tol: + keep = False + break + if keep: + indices_to_keep.add(i) + filtered_peaks = peaks[list(indices_to_keep)] + + vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp))) # this tapers the start + nmax = np.argmax(ht.data.data) ht.data.data[0:n_samp] *= vectaper - # If Scattering waveform, tapers same amount as early taper - # If Capture waveform, no end taper - if np.abs(ht.data.length-nmax) > 3e3: + + if len(filtered_peaks) == 1: + # check if a scatter or a plunge + if np.abs(ht.data.length-nmax) > 3e3: + #scatter + print('Tapering for scatter waveform') + n_samp2=n_samp + vectaper2 = 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2)) # end taper + ht.data.data[-(n_samp2+1):] *= vectaper2 + else: + # plunge, only need to taper the start + print('Plunge waveform, only start taper') + + elif len(filtered_peaks) == 0: + print('Non-interactive hyperbolic waveform, tapering both ends') n_samp2=n_samp - vectaper2= 0.5 - 0.5*np.cos(np.pi* (1-np.arange(n_samp2+1)/(1.*n_samp2))) + vectaper2 = 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2)) ht.data.data[-(n_samp2+1):] *= vectaper2 + + else: + # zoom-whirl case, taper just the start + print('Zoom-whirl waveform, only start taper') + if P.deltaF is not None: TDlen = int(1./P.deltaF * 1./P.deltaT) @@ -3029,7 +3307,7 @@ def hoft(P, Fp=None, Fc=None,**kwargs): ht = lal.ResizeREAL8TimeSeries(ht, 0, TDlen) return ht -def hoff(P, Fp=None, Fc=None, fwdplan=None): +def hoff(P, Fp=None, Fc=None, fwdplan=None, **kwargs): """ Generate a FD waveform from ChooseWaveformParams P. Will return a COMPLEX16FrequencySeries object. @@ -3045,6 +3323,8 @@ def hoff(P, Fp=None, Fc=None, fwdplan=None): If P.deltaF == None, the TD waveform will be zero-padded to the next power of 2. """ + Lmax = kwargs.get('Lmax', 4) # Default to 4 if not specified in kwargs + # For FD approximants, use the ChooseFDWaveform path = hoff_FD if lalsim.SimInspiralImplementedFDApproximants(P.approx)==1: # Raise exception if unused arguments were specified @@ -3054,11 +3334,11 @@ def hoff(P, Fp=None, Fc=None, fwdplan=None): # For TD approximants, do ChooseTDWaveform + FFT path = hoff_TD else: - hf = hoff_TD(P, Fp, Fc, fwdplan) + hf = hoff_TD(P, Fp, Fc, fwdplan, Lmax=Lmax) return hf -def hoff_TD(P, Fp=None, Fc=None, fwdplan=None): +def hoff_TD(P, Fp=None, Fc=None, fwdplan=None, **kwargs): """ Generate a FD waveform from ChooseWaveformParams P by creating a TD waveform, zero-padding and @@ -3078,7 +3358,9 @@ def hoff_TD(P, Fp=None, Fc=None, fwdplan=None): Returns a COMPLEX16FrequencySeries object """ - ht = hoft(P, Fp, Fc) + Lmax = kwargs.get('Lmax', 4) # Default to 4 if not specified in kwargs + + ht = hoft(P, Fp, Fc, Lmax=Lmax) if P.deltaF == None: # h(t) was not zero-padded, so do it now TDlen = nextPow2(ht.data.length) @@ -3438,6 +3720,7 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil M1=P.m1/lal.MSUN_SI M2=P.m2/lal.MSUN_SI nu=M1*M2/((M1+M2)**2) + hyp_wav = False if (P.eccentricity == 0.0 and P.E0 ==0.0): print("Using ResumS master; not eccentric") pars = { @@ -3467,6 +3750,10 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil elif (P.eccentricity == 0.0): print("Using hyperbolic call RIFT O4b branch") + hyp_wav = True # convenient way to know if the waveform is hyperbolic + if kwargs.get('force_22_mode', False): + print('Forcing ONLY the 22 modes') + k = [1] pars = { 'M' : M1+M2, 'q' : M1/M2, @@ -3527,7 +3814,102 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil t, hptmp, hctmp, hlmtmp, dym = EOBRun_module.EOBRunPy(pars) print("EOBRun_module done") k_list_orig = hlmtmp.keys() - hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2) + + if not(hyp_wav): + ## Set the epoch for non-hyperbolic cases ## + hpepoch = -P.deltaT*np.argmax(np.abs(hptmp)**2+np.abs(hctmp)**2) + else: + ## custom epoch for the hyperbolic case ## + + # wf amplitude - 22 mode only + tmp_22 = np.array(hlmtmp['1']) + tmp_22[0] *= (m_total_s/distance_s)*nu + amp = np.abs(tmp_22[0] * np.exp(-1j*(2*(np.pi/2.)+tmp_22[1]))) + amp_norm = amp / np.amax(amp) # normalize amplitude for peak finding + + # Check for cases where the amplitude is max at the start or end + if np.argmax(amp_norm) == 0 or np.argmax(amp_norm) == len(amp_norm) - 1: + reclassify = True + else: + reclassify = False + + # initial peak finding to determine system type + height_thresh = 0.25 + prom_thresh = 0.1 + peaks, props = signal.find_peaks(amp_norm, height = height_thresh, prominence = prom_thresh) + peak_heights = props['peak_heights'] + # filtering out peaks so we only keep the local maxima + indices_to_keep = set() + sorted_indices = np.argsort(peak_heights)[::-1] + tol = int(pars['srate_interp'] / 13.65) # 300 samples at srate of 4096 - minimum distance between peaks. + for i in sorted_indices: + peak = peaks[i] + keep = True + for kept_index in indices_to_keep: + if abs(peaks[kept_index] - peak) <= tol: + keep = False + break + if keep: + indices_to_keep.add(i) + filtered_peaks = peaks[list(indices_to_keep)] + # parsing number of peaks after filtering against distance tolerance + if len(filtered_peaks) == 1: + # scatter case OR plunge case, we can set the epoch normally + hpepoch = -P.deltaT*np.argmax(amp) + + if np.abs(amp)[-1] > 1e-26: + hypclass = 'scatter' # maybe should do through assign_param? + else: + hypclass = 'plunge' + + + elif len(filtered_peaks) == 0: + hypclass = 'meaningless' + hpepoch = -P.deltaT*np.argmax(amp) + reclassify = True + + else: + # capture case, we need to force the epoch to be the last peak + hypclass = 'zoomwhirl' + hpepoch = -P.deltaT*filtered_peaks[-1] + + if reclassify: + print('Running re-classification') + # run minimal threshold peak finder + all_peaks, all_props = signal.find_peaks(amp_norm, height=0.000001, prominence=0.000001) + + # checking for minima if no peaks found + + if len(all_peaks) == 0: + print("No peaks found, checking for minima instead...") + all_peaks, all_props = signal.find_peaks((-1*amp_norm + 1.0), height=0.000001, prominence=0.000001) + + if len(all_props['prominences']) > 3: + print('MANY peaks detected on reclassification, evaluating...') + + if np.abs(amp)[-1] < 1e-26: + print('Reclassifying to Plunge') + hpepoch = -P.deltaT*np.argmax(amp) + hypclass = 'plunge' + else: + print('Reclassifying to Scatter') + max_peak_index = all_peaks[np.argmax(all_props['peak_heights'])] + print(f"Largest peak is at index {max_peak_index} with height {amp_norm[max_peak_index]}") + # setting the epoch to the largest amplitude peak + hpepoch = -P.deltaT*max_peak_index + elif len(all_props['prominences']) == 3 or len(all_props['prominences']) == 2 or len(all_props['prominences']) == 1: + print('Reclassifying to Scatter') + hypclass = 'scatter' + + max_peak_index = all_peaks[np.argmax(all_props['peak_heights'])] + print(f"Largest peak is at index {max_peak_index} with height {amp_norm[max_peak_index]}") + + # setting the epoch to the largest amplitude peak + hpepoch = -P.deltaT*max_peak_index + else: + print('No peaks detected after reclassifcation') + hypclass='meaningless' + hlmlen = len(hptmp) hlm = {} hlmtmp2 = {} @@ -3599,7 +3981,8 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil modes_used_new.append((4,0)) hlmtmp2[(4,0)]=np.array(hlmtmp[k]) modes_used=modes_used_new - print(modes_used,hlmtmp,hlmtmp2) + if not(hyp_wav): + print(modes_used,hlmtmp,hlmtmp2) # for count,mode in enumerate(modes_used): # hlmtmp2[mode]=np.array(hlmtmp[str(count)]) for mode in modes_used: @@ -3639,6 +4022,7 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil for mode in modes_used_new2: hlm[mode].data.data[0:n_samp] *= vectaper else: + # determine location of start taper if not 'n_samp' in locals(): for count,value in enumerate(hlm[(2,2)].data.data): if count ==0: @@ -3646,17 +4030,42 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil if np.abs(np.real(value)-np.real(hlm[(2,2)].data.data[0])) > 0.01 * np.abs(np.real(hlm[(2,2)].data.data[0])): n_samp=int(count/2) break + + # determine location of end taper + if not 'n_samp2' in locals(): + for count, value in enumerate(reversed(hlm[(2,2)].data.data)): # Scan backwards + if count == 0: + continue + if np.abs(np.real(value) - np.real(hlm[(2,2)].data.data[-1])) > 0.01 * np.abs(np.real(hlm[(2,2)].data.data[-1])): + n_samp2 = int(count / 2) + break + + # Always taper the start + vectaper= 0.5 + 0.5*np.cos(np.pi* (1-np.arange(n_samp)/(1.*n_samp))) nmax = np.argmax(hlm[(2,2)].data.data) for mode in modes_used_new2: + #pass hlm[mode].data.data[0:n_samp] *= vectaper - # If Capture waveform, no end taper - # If Scattering waveform, tapers end same amount as early taper - if np.abs(hlm[(2,2)].data.length-nmax) > 3e3: - n_samp2=n_samp - vectaper2= 0.5 - 0.5*np.cos(np.pi* (1-np.arange(n_samp2+1)/(1.*n_samp2))) + + if hypclass == 'scatter': + # Taper for scatter + vectaper2= 0.5 + 0.5 * np.cos(np.pi * np.arange(n_samp2 + 1) / (1. * n_samp2)) for mode in modes_used_new2: hlm[mode].data.data[-(n_samp2+1):] *= vectaper2 + elif hypclass == 'plunge': + # taper for plunge + print('Plunge waveform, only start taper') + elif hypclass == 'zoomwhirl': + # taper for ZW + print('Zoom-whirl waveform, only start taper') + elif hypclass =='meaningless': + # zero out meaningless + for mode in modes_used_new2: + hlm[mode].data.data *= 0.0 + + # end off custom hyperbolic stuff + for mode in modes_used_new2: if not (P.deltaF is None): TDlen = int(1./P.deltaF * 1./P.deltaT) @@ -3676,6 +4085,9 @@ def hlmoft(P, Lmax=2,nr_polarization_convention=False, fixed_tapering=False, sil # for mode in hlm: # hlm[mode].data.data*= hp.data.data + print('Actual modes used:') + print(hlm.keys()) + return hlm else: # (P.approx == lalSEOBv4 or P.approx == lalsim.SEOBNRv2 or P.approx == lalsim.SEOBNRv1 or P.approx == lalsim.EOBNRv2 extra_params = P.to_lal_dict_extended(extra_args_dict=extra_waveform_args) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py index 6b3b1222..df25912b 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py @@ -121,7 +121,7 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict, extra_waveform_kwargs={}, use_gwsignal=False, use_gwsignal_approx=None, - use_external_EOB=False,nr_lookup=False,nr_lookup_valid_groups=None,no_memory=True,perturbative_extraction=False,perturbative_extraction_full=False,hybrid_use=False,hybrid_method='taper_add',use_provided_strain=False,ROM_group=None,ROM_param=None,ROM_use_basis=False,ROM_limit_basis_size=None,skip_interpolation=False): + use_external_EOB=False,nr_lookup=False,nr_lookup_valid_groups=None,no_memory=True,perturbative_extraction=False,perturbative_extraction_full=False,hybrid_use=False,hybrid_method='taper_add',use_provided_strain=False,ROM_group=None,ROM_param=None,ROM_use_basis=False,ROM_limit_basis_size=None,skip_interpolation=False,**kwargs): """ Compute < h_lm(t) | d > and < h_lm | h_l'm' > @@ -142,6 +142,10 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict, rholms_intp = {} crossTerms = {} crossTermsV = {} + + if kwargs.get('force_22_mode', False): + # forcing 22 mode only for hyperbolic case + force_22_mode = True # Compute hlms at a reference distance, distance scaling is applied later P.dist = distMpcRef*1e6*lsu.lsu_PC @@ -198,7 +202,7 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict, elif use_gwsignal and (has_GWS): # this MUST be called first, so the P.approx is never tested if not quiet: print( " FACTORED LIKELIHOOD WITH hlmoff (GWsignal) " ) - hlms, hlms_conj = rgws.std_and_conj_hlmoff(P,Lmax,approx_string=use_gwsignal_approx,**extra_waveform_kwargs) + hlms, hlms_conj = rgws.std_and_conj_hlmoff(P,Lmax,approx_string=use_gwsignal_approx,force_22_mode=force_22_mode,**extra_waveform_kwargs) elif (not nr_lookup) and (not NR_group) and ( P.approx ==lalsim.SEOBNRv2 or P.approx == lalsim.SEOBNRv1 or P.approx==lalsim.SEOBNRv3 or P.approx == lsu.lalSEOBv4 or P.approx ==lsu.lalSEOBNRv4HM or P.approx == lalsim.EOBNRv2 or P.approx == lsu.lalTEOBv2 or P.approx==lsu.lalTEOBv4 ): # note: alternative to this branch is to call hlmoff, which will actually *work* if ChooseTDModes is propertly implemented for that model @@ -250,7 +254,7 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict, # hlms_conj = lsu.SphHarmFrequencySeries_to_dict(hlms_conj_list, Lmax) # a dictionary # else: # hlms_conj = hlms_conj_list - hlms, hlms_conj = lsu.std_and_conj_hlmoff(P,Lmax,**extra_waveform_kwargs) + hlms, hlms_conj = lsu.std_and_conj_hlmoff(P,Lmax,force_22_mode=force_22_mode,**extra_waveform_kwargs) elif (nr_lookup or NR_group) and useNR: # look up simulation # use nrwf to get hlmf diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/samples_utils.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/samples_utils.py index 4d96e57f..06428531 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/misc/samples_utils.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/samples_utils.py @@ -137,6 +137,18 @@ def extract_combination_from_LI(samples_LI, p): a = a[8:] if a =='q' and 'm1' in samples_LI.dtype.names: return samples_LI["m1"]/samples_LI["m2"] + + if p == 'b_hyp': + samples = np.array([samples_LI["m1"], samples_LI["m2"], samples_LI["a1x"], samples_LI["a1y"], samples_LI["a1z"], samples_LI["a2x"], samples_LI["a2y"], samples_LI["a2z"], samples_LI["E0"], samples_LI["p_phi0"]]).T + with Pool(12) as pool: + b_hyp = np.array(pool.map(f_b_hyp, samples)) + return b_hyp + + if p == 'phi_scatter': + samples = np.array([samples_LI["m1"], samples_LI["m2"], samples_LI["a1x"], samples_LI["a1y"], samples_LI["a1z"], samples_LI["a2x"], samples_LI["a2y"], samples_LI["a2z"], samples_LI["E0"], samples_LI["p_phi0"]]).T + with Pool(12) as pool: + phi_scatter = np.array(pool.map(f_phi_scatter, samples)) + return phi_scatter print(" No access for parameter ", p) return np.zeros(len(samples_LI['m1'])) # to avoid causing a hard failure @@ -272,3 +284,33 @@ def fchip(sample): P.s2z = sample[7] chip = P.extract_param('chi_p') return chip + +def f_b_hyp(sample): + P=lalsimutils.ChooseWaveformParams() + P.m1 = sample[0] + P.m2 = sample[1] + P.s1x = sample[2] + P.s1y = sample[3] + P.s1z = sample[4] + P.s2x = sample[5] + P.s2y = sample[6] + P.s2z = sample[7] + P.E0 = sample[8] + P.p_phi0 = sample[9] + b_hyp = P.extract_param('b_hyp') + return b_hyp + +def f_phi_scatter(sample): + P=lalsimutils.ChooseWaveformParams() + P.m1 = sample[0] + P.m2 = sample[1] + P.s1x = sample[2] + P.s1y = sample[3] + P.s1z = sample[4] + P.s2x = sample[5] + P.s2y = sample[6] + P.s2z = sample[7] + P.E0 = sample[8] + P.p_phi0 = sample[9] + phi_scatter = P.extract_param('phi_scatter') + return phi_scatter \ No newline at end of file diff --git a/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py b/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py index 54fffc65..c6ce3976 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/misc/xmlutils.py @@ -42,6 +42,8 @@ def assign_time(row, t): "joint_prior": "alpha2", "joint_s_prior": "alpha3", "eccentricity":"alpha4", + "E0":"psi3", + "p_phi0":"beta", "lambda1":"alpha5", "lambda2":"alpha6", "spin1x":"spin1x", diff --git a/MonteCarloMarginalizeCode/Code/bin/convert_output_format_ile2inference b/MonteCarloMarginalizeCode/Code/bin/convert_output_format_ile2inference index 6927cf43..8ef5d4a3 100755 --- a/MonteCarloMarginalizeCode/Code/bin/convert_output_format_ile2inference +++ b/MonteCarloMarginalizeCode/Code/bin/convert_output_format_ile2inference @@ -39,7 +39,7 @@ optp.add_option("--export-eos",action='store_true',help="Include EOS index param optp.add_option("--export-cosmology",action='store_true',help="Include source frame masses and redshift") optp.add_option("--export-weights",action='store_true',help="Include a field 'weights' equal to L p/ps") optp.add_option("--export-eccentricity", action="store_true", help="Include eccentricity") -optp.add_option("--export-hyperbolic", action="store_true", help="Include eccentricity") +optp.add_option("--export-hyperbolic", action="store_true", help="Include E0 and p_phi0") optp.add_option("--with-cosmology",default="Planck15",help="Specific cosmology to use") optp.add_option("--use-interpolated-cosmology",action='store_true',help="Specific cosmology to use") optp.add_option("--convention",default="RIFT",help="RIFT|LI") @@ -99,10 +99,11 @@ if opts.convention == 'LI': if opts.export_eccentricity: ecc = [row.alpha4 for row in points] + if opts.export_hyperbolic: p_phi0 = [row.beta for row in points] E0 = [row.psi3 for row in points] - + wt = np.exp(like)*p/ps neff_here = np.sum(wt)/np.max(wt) # neff for this file. Assumes NOT mixed samples: dangerous diff --git a/MonteCarloMarginalizeCode/Code/bin/convert_output_format_inference2ile b/MonteCarloMarginalizeCode/Code/bin/convert_output_format_inference2ile index f84696cb..077fe590 100755 --- a/MonteCarloMarginalizeCode/Code/bin/convert_output_format_inference2ile +++ b/MonteCarloMarginalizeCode/Code/bin/convert_output_format_inference2ile @@ -84,6 +84,10 @@ for indx in np.arange(opts.target_size): P = lalsimutils.ChooseWaveformParams(m1=m1,m2=m2,dist=d) if "eccentricity" in samples_in.dtype.names: P.eccentricity = samples_in["eccentricity"][fac_reduce*indx] + if "E0" in samples_in.dtype.names: + P.E0 = samples_in["E0"][fac_reduce*indx] + if "p_phi0" in samples_in.dtype.names: + P.p_phi0 = samples_in["p_phi0"][fac_reduce*indx] if "time_maxl" in samples_in.dtype.names: P.tref = samples_in["time_maxl"][fac_reduce*indx] elif 'time' in samples_in.dtype.names: diff --git a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py index 707cad87..36ed9bd5 100755 --- a/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py +++ b/MonteCarloMarginalizeCode/Code/bin/helper_LDG_Events.py @@ -157,6 +157,7 @@ def get_observing_run(t): parser.add_argument("--limit-mc-range",default=None,type=str,help="For PP plots, or other analyses requiring a specific mc range (eg ini file), bounding the limit *above*. Allows the code to auto-select its mc range as usual, then takes the intersection with this limit") parser.add_argument("--scale-mc-range",type=float,default=None,help="If using the auto-selected mc, scale the ms range proposed by a constant factor. Recommend > 1. . ini file assignment will override this.") parser.add_argument("--force-eta-range",default=None,type=str,help="For PP plots. Enforces initial grid placement inside this region") +parser.add_argument("--force-mtot-range",default=None,type=str,help="For PP plots, hyperbolic analysis, or other analyses requiring a specific mtot range (eg ini file). Enforces initial grid placement inside this region. Passed directly to MOG and CIP.") parser.add_argument("--allow-subsolar", action='store_true', help="Override limits which otherwise prevent subsolar mass PE") parser.add_argument("--use-legacy-gracedb",action='store_true') parser.add_argument("--event-time",type=float,default=None) @@ -198,7 +199,12 @@ def get_observing_run(t): parser.add_argument("--assume-matter-but-primary-bh",action='store_true',help="If present, the code will add options necessary to manage tidal arguments for the smaller body ONLY. (Usually pointless)") parser.add_argument("--internal-tabular-eos-file",type=str,default=None,help="Tabular file of EOS to use. The default prior will be UNIFORM in this table!. NOT YET IMPLEMENTED (initial grids, etc)") parser.add_argument("--assume-eccentric",action='store_true',help="If present, the code will add options necessary to manage eccentric arguments. The proposed fit strategy and initial grid will allow for eccentricity") -parser.add_argument("--assume-hyperbolic",action='store_true',help="If present, the code will add options necessary to manage eccentric arguments. The proposed fit strategy and initial grid will allow for hyperbolic") +parser.add_argument("--assume-hyperbolic",action='store_true',help="If present, the code will add options necessary to manage hyperbolic arguments. The proposed fit strategy and initial grid will allow for hyperbolic orbits") +parser.add_argument("--E0-max", default=1.060,type=float,help="Maximum range of 'E0' allowed.") +parser.add_argument("--E0-min", default=1.0,type=float,help="Minimum range of 'E0' allowed.") +parser.add_argument("--pphi0-max", default=5.4,type=float,help="Maximum range of 'p_phi0' allowed.") +parser.add_argument("--pphi0-min", default=3.8,type=float,help="Minumum range of 'p_phi0' allowed.") +parser.add_argument("--use-mtot-coords",action='store_true',help="Configures CIP and PUFF for mtot instead of mc. REQUIRES --force-mtot-range.") parser.add_argument("--assume-nospin",action='store_true',help="If present, the code will not add options to manage precessing spins (the default is aligned spin)") parser.add_argument("--assume-precessing-spin",action='store_true',help="If present, the code will add options to manage precessing spins (the default is aligned spin)") parser.add_argument("--assume-volumetric-spin",action='store_true',help="If present, the code will assume a volumetric spin prior in its last iterations. If *not* present, the code will adopt a uniform magnitude spin prior in its last iterations. If not present, generally more iterations are taken.") @@ -246,8 +252,27 @@ def get_observing_run(t): parser.add_argument("--use-cvmfs-frames",action='store_true',help="If true, require LIGO frames are present (usually via CVMFS). User is responsible for generating cache file compatible with it. This option insures that the cache file is properly transferred (because you have generated it)") parser.add_argument("--use-ini",default=None,type=str,help="Attempt to parse LI ini file to set corresponding options. WARNING: MAY OVERRIDE SOME OTHER COMMAND-LINE OPTIONS") parser.add_argument("--verbose",action='store_true') +parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.") +parser.add_argument("--force-plunge-grids",action='store_true',help="Eliminates all non-plunge intrinsic points from hyperbolic grids throughout the workflow.") +parser.add_argument("--force-zoomwhirl-grids",action='store_true',help="Eliminates all non-zoomwhirl intrinsic points from hyperbolic grids throughout the workflow.") +parser.add_argument('--force-hyperbolic-22', action='store_true', help='Forces just the 22 modes for hyperbolic waveforms') opts= parser.parse_args() +# Ensure --assume-hyperbolic is set when using any --force-X-grids option +# Ensure only ONE of the --force-X-grids options is set +force_grids = [opts.force_scatter_grids, opts.force_plunge_grids, opts.force_zoomwhirl_grids] +if any(force_grids) and not opts.assume_hyperbolic: + parser.error("Using --force-scatter-grids, --force-plunge-grids, or --force-zoomwhirl-grids requires --assume-hyperbolic!") + +if sum(bool(x) for x in force_grids) > 1: + parser.error("CANNOT use multiple --force-X-grids options at the same time!") + +if opts.use_mtot_coords: + if opts.force_mtot_range is None: + print('Using the mtot coords requires a specified range!') + print('Specify the mtot range with --force-mtot-range!') + sys.exit(1) + if opts.assume_matter_but_primary_bh: opts.assume_matter=True @@ -451,6 +476,9 @@ def get_observing_run(t): event_dict["s2z"] = P.s2z event_dict["P"] = P event_dict["epoch"] = 0 # no estimate for now + if opts.assume_hyperbolic: + event_dict["E0"] = P.E0 + event_dict["p_phi0"] = P.p_phi0 elif opts.use_coinc: # If using a coinc through injections and not a GraceDB event. # Same code as used before for gracedb coinc_file = opts.use_coinc @@ -980,10 +1008,6 @@ def crit_m2(delta): eta_min = q_min/(1.+q_min)**2 if 'ecc_min' in engine_dict: ecc_range_str = " ["+str(engine_dict['force_ecc_min'])+","+str(engine_dict['force_ecc_max'])+"]" - if 'E0_min' in engine_dict: - E0_range_str = " ["+str(engine_dict['force_E0_min'])+","+str(engine_dict['force_E0_max'])+"]" - if 'pphi0_min' in engine_dict: - pphi0_range_str = " ["+str(engine_dict['force_pphi0_min'])+","+str(engine_dict['force_pphi0_max'])+"]" mc_range_str = " ["+str(mc_min_tight)+","+str(mc_max_tight)+"]" # Use a tight placement grid for CIP if not(opts.manual_mc_min is None): @@ -993,6 +1017,9 @@ def crit_m2(delta): mc_range_str_cip = " --mc-range ["+str(mc_min)+","+str(mc_max)+"]" if not(opts.force_mc_range is None): mc_range_str_cip = " --mc-range " + opts.force_mc_range +if not(opts.force_mtot_range is None): + mtot_range_str = opts.force_mtot_range + mtot_range_str_cip = " --mtot-range " + opts.force_mtot_range elif opts.limit_mc_range: line_here = list(map(float,opts.limit_mc_range.replace('[','').replace(']','').split(',') )) mc_min_lim,mc_max_lim = line_here @@ -1156,16 +1183,24 @@ def crit_m2(delta): helper_ile_args += " --data-start-time " + str(data_start_time) + " --data-end-time " + str(data_end_time) + " --inv-spec-trunc-time 0 --window-shape " + str(window_shape) elif opts.data_LI_seglen: seglen = opts.data_LI_seglen - # Use LI-style positioning of trigger relative to 2s before end of buffer - # Use LI-style tukey windowing - window_shape = 0.4*2/seglen - data_end_time = event_dict["tref"]+2 - data_start_time = event_dict["tref"] +2 - seglen + if opts.assume_hyperbolic: + window_shape = 0.0 # DO NOT window at all + # split seglen across the event time + data_start_time = event_dict["tref"] - seglen/2 + data_end_time = event_dict["tref"] + seglen/2 + else: + # Use LI-style positioning of trigger relative to 2s before end of buffer + # Use LI-style tukey windowing + window_shape = 0.4*2/seglen + data_end_time = event_dict["tref"]+2 + data_start_time = event_dict["tref"] +2 - seglen helper_ile_args += " --data-start-time " + str(data_start_time) + " --data-end-time " + str(data_end_time) + " --inv-spec-trunc-time 0 --window-shape " + str(window_shape) if opts.assume_eccentric: helper_ile_args += " --save-eccentricity " if opts.assume_hyperbolic: helper_ile_args += " --save-hyperbolic " + if opts.force_hyperbolic_22: + helper_ile_args += " --force-hyperbolic-22 " if opts.propose_initial_grid_fisher: # and (P.extract_param('mc')/lal.MSUN_SI < 10.): cmd = "util_AnalyticFisherGrid.py --inj-file-out proposed-grid " # Add standard downselects : do not have m1, m2 be less than 1 @@ -1196,8 +1231,7 @@ def crit_m2(delta): if opts.assume_eccentric: cmd += " --random-parameter eccentricity --random-parameter-range " + ecc_range_str if opts.assume_hyperbolic: - cmd += " --random-parameter E0 --random-parameter-range " + E0_range_str - cmd += " --random-parameter p_phi0 --random-parameter-range " + pphi0_range_str + cmd += f" --random-parameter E0 --random-parameter-range [{opts.E0_min},{opts.E0_max}] --random-parameter p_phi0 --random-parameter-range [{opts.pphi0_min},{opts.pphi0_max}] " if "SNR" in event_dict: grid_size *= np.max([1,event_dict["SNR"]/15]) # more grid points at higher amplitude. Yes, even though we also contract the paramete range if not (opts.force_initial_grid_size is None): @@ -1214,13 +1248,18 @@ def crit_m2(delta): # add basic mass parameters cmd = "util_ManualOverlapGrid.py --fname proposed-grid --skip-overlap " mass_string_init = " --random-parameter mc --random-parameter-range " + mc_range_str + " --random-parameter delta_mc --random-parameter-range '[" + str(delta_grid_min) +"," + str(delta_grid_max) + "]' " + if not(opts.force_mtot_range is None): + mass_string_init = " --random-parameter mtot --random-parameter-range " + mtot_range_str + " --random-parameter delta_mc --random-parameter-range '[" + str(delta_grid_min) +"," + str(delta_grid_max) + "]' " cmd+= mass_string_init # Add standard downselects : do not have m1, m2 be less than 1 if not(opts.force_mc_range is None): # force downselect based on this range cmd += " --downselect-parameter mc --downselect-parameter-range " + opts.force_mc_range if not(opts.force_eta_range is None): - cmd += " --downselect-parameter eta --downselect-parameter-range " + opts.force_eta_range + cmd += " --downselect-parameter eta --downselect-parameter-range " + opts.force_eta_range + if not(opts.force_mtot_range is None): + # force downselect based on this range + cmd += " --downselect-parameter mtot --downselect-parameter-range " + opts.force_mtot_range cmd += " --fmin " + str(opts.fmin_template) if opts.data_LI_seglen and not (opts.no_enforce_duration_bound): cmd += " --enforce-duration-bound " + str(opts.data_LI_seglen) @@ -1253,8 +1292,13 @@ def crit_m2(delta): if opts.assume_eccentric: cmd += " --random-parameter eccentricity --random-parameter-range " + ecc_range_str if opts.assume_hyperbolic: - cmd += " --random-parameter E0 --random-parameter-range " + E0_range_str - cmd += " --random-parameter p_phi0 --random-parameter-range " + pphi0_range_str + cmd += f" --random-parameter E0 --random-parameter-range [{opts.E0_min},{opts.E0_max}] --random-parameter p_phi0 --random-parameter-range [{opts.pphi0_min},{opts.pphi0_max}] " + if opts.force_scatter_grids: + cmd += " --force-scatter " + if opts.force_plunge_grids: + cmd += " --force-plunge " + if opts.force_zoomwhirl_grids: + cmd += " --force-zoomwhirl " if opts.internal_tabular_eos_file: cmd += " --tabular-eos-file {} ".format(opts.internal_tabular_eos_file) grid_size *=2 # larger grids needed for discrete realization scenarios @@ -1304,6 +1348,7 @@ def lambda_m_estimate(m): if not (opts.force_initial_grid_size is None): grid_size = opts.force_initial_grid_size + cmd += " --grid-cartesian-npts " + str(int(grid_size)) print(" Executing grid command ", cmd) os.system(cmd) @@ -1398,14 +1443,24 @@ def lambda_m_estimate(m): puff_max_it=0 helper_puff_args = " --parameter mc --parameter eta --fmin {} --fref {} ".format(opts.fmin_template,opts.fmin_template) +if opts.use_mtot_coords: + helper_puff_args = " --parameter mtot --parameter q --fmin {} --fref {} ".format(opts.fmin_template,opts.fmin_template) if opts.assume_eccentric: helper_puff_args += " --parameter eccentricity " if opts.assume_hyperbolic: helper_puff_args += " --parameter E0 --parameter p_phi0 " + if opts.force_scatter_grids: + helper_puff_args += " --force-scatter " + if opts.force_plunge_grids: + helper_puff_args += " --force-plunge " + if opts.force_zoomwhirl_grids: + helper_puff_args += " --force-zoomwhirl " if event_dict["MChirp"] >25: # at high mass, mc/eta correlation weak, don't want to have eta coordinate degeneracy at q=1 to reduce puff proposals near there - helper_puff_args = " --parameter mc --parameter delta_mc " + helper_puff_args = " --parameter mc --parameter delta_mc " + if opts.use_mtot_coords: + helper_puff_args = " --parameter mtot --parameter delta_mc " if opts.propose_fit_strategy: puff_max_it= 0 # Strategy: One iteration of low-dimensional, followed by other dimensions of high-dimensional @@ -1419,7 +1474,10 @@ def lambda_m_estimate(m): if 'gp' in fit_method: helper_cip_args += " --cap-points 12000 " if not opts.no_propose_limits: - helper_cip_args += mc_range_str_cip + eta_range_str_cip + if not(opts.use_mtot_coords): + helper_cip_args += mc_range_str_cip + eta_range_str_cip + else: + helper_cip_args += mtot_range_str_cip + eta_range_str_cip if opts.force_chi_max: helper_cip_args += " --chi-max {} ".format(opts.force_chi_max) if opts.force_chi_small_max: @@ -1722,6 +1780,7 @@ def lambda_m_estimate(m): if opts.assume_eccentric: with open("helper_convert_args.txt",'w+') as f: f.write(" --export-eccentricity ") + if opts.assume_hyperbolic: with open("helper_convert_args.txt",'w+') as f: f.write(" --export-hyperbolic ") diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode index c132cc37..9829927e 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode @@ -241,6 +241,7 @@ optp.add_option("--internal-waveform-extra-lalsuite-args",type=str,default=None) optp.add_option("--verbose",action='store_true') optp.add_option("--save-eccentricity", action="store_true") optp.add_option("--save-hyperbolic", action="store_true") +optp.add_option('--force-hyperbolic-22', default=False, action='store_true', help='Forces just the 22 modes for hyperbolic waveforms') # # Add the integration options # @@ -1187,7 +1188,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t NR_group=NR_template_group,NR_param=NR_template_param, use_gwsignal=opts.use_gwsignal, use_gwsignal_approx=opts.approximant, - use_external_EOB=opts.use_external_EOB,nr_lookup=opts.nr_lookup,nr_lookup_valid_groups=opts.nr_lookup_group,perturbative_extraction=opts.nr_perturbative_extraction,perturbative_extraction_full=opts.nr_perturbative_extraction_full,use_provided_strain=opts.nr_use_provided_strain,hybrid_use=opts.nr_hybrid_use,hybrid_method=opts.nr_hybrid_method,ROM_group=opts.rom_group,ROM_param=opts.rom_param,ROM_use_basis=opts.rom_use_basis,verbose=opts.verbose,quiet=not opts.verbose,ROM_limit_basis_size=opts.rom_limit_basis_size_to,no_memory=opts.no_memory,skip_interpolation=opts.vectorized, extra_waveform_kwargs=extra_waveform_kwargs) + use_external_EOB=opts.use_external_EOB,nr_lookup=opts.nr_lookup,nr_lookup_valid_groups=opts.nr_lookup_group,perturbative_extraction=opts.nr_perturbative_extraction,perturbative_extraction_full=opts.nr_perturbative_extraction_full,use_provided_strain=opts.nr_use_provided_strain,hybrid_use=opts.nr_hybrid_use,hybrid_method=opts.nr_hybrid_method,ROM_group=opts.rom_group,ROM_param=opts.rom_param,ROM_use_basis=opts.rom_use_basis,verbose=opts.verbose,quiet=not opts.verbose,ROM_limit_basis_size=opts.rom_limit_basis_size_to,no_memory=opts.no_memory,skip_interpolation=opts.vectorized, extra_waveform_kwargs=extra_waveform_kwargs, force_22_mode=opts.force_hyperbolic_22) if opts.auto_logarithm_offset and guess_snr: # important: this only impacts *this* analysis @@ -1636,7 +1637,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, P.eccentricity, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]])) #dict_return["convergence_test_results"]["normal_integral]" elif opts.save_hyperbolic: # output format when hyperbolic is being used - numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, P.E0, P.p_phi0, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]])) #dict_return["convergence_test_results"]["normal_integral]" + numpy.savetxt(fname_output_txt, numpy.array([[event_id, m1, m2, P.s1x, P.s1y, P.s1z, P.s2x, P.s2y, P.s2z, P.E0, P.p_phi0, log_res+manual_avoid_overflow_logarithm, sqrt_var_over_res,sampler.ntotal, neff ]])) elif not (P.lambda1>0 or P.lambda2>0): # output format when lambda is NOT used if not opts.pin_distance_to_sim: @@ -1716,7 +1717,6 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t samples["alpha6"] =numpy.ones(samples["psi"].shape)*P.lambda2 samples["psi3"] =numpy.ones(samples["psi"].shape)*P.E0 samples["beta"] =numpy.ones(samples["psi"].shape)*P.p_phi0 - # Below exist solely to placate XML export; new issue as of latest lalsuite say 7.15+ or so samples["alpha2"] = numpy.zeros(samples["psi"].shape) samples["alpha3"] = numpy.zeros(samples["psi"].shape) @@ -1727,7 +1727,7 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t # print(samples['t_ref'] - fiducial_epoch, len(samples['t_ref'])) xmlutils.append_samples_to_xmldoc(xmldoc, samples) # Extra metadata - dict_out={"mass1": m1, "mass2": m2, "spin1z": P.s1z, "spin2z": P.s2z, "alpha4": P.eccentricity, "alpha5":P.lambda1, "alpha6":P.lambda2, "psi3":P.E0, "beta":P.p_phi0, "event_duration": sqrt_var_over_res, "ttotal": sampler.ntotal} + dict_out={"mass1": m1, "mass2": m2, "spin1z": P.s1z, "spin2z": P.s2z, "alpha4": P.eccentricity, "alpha5":P.lambda1, "alpha6":P.lambda2, "psi3":P.E0, "beta":P.p_phi0, "event_duration": sqrt_var_over_res, "ttotal": sampler.ntotal} # if 'distance' in pinned_params: # dict_out['distance'] = pinned_params["distance"] converged_result = False diff --git a/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py b/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py index a1fd60d3..91ed35b2 100755 --- a/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py +++ b/MonteCarloMarginalizeCode/Code/bin/plot_posterior_corner.py @@ -232,6 +232,7 @@ def render_coordinates(coord_names,logparams=[]): parser.add_argument("--lnL-cut",default=None,type=float) parser.add_argument("--sigma-cut",default=0.4,type=float) parser.add_argument("--eccentricity", action="store_true", help="Read sample files in format including eccentricity") +parser.add_argument("--hyperbolic", action="store_true", help="Read sample files in format including hyperbolic params") parser.add_argument("--matplotlib-block-defaults",action="store_true",help="Relies entirely on user to set plot options for plot styles from matplotlibrc") parser.add_argument("--no-mod-psi",action="store_true",help="Default is to take psi mod pi. If present, does not do this") parser.add_argument("--verbose",action='store_true',help='print matplotlibrc data') @@ -467,6 +468,9 @@ def render_coordinates(coord_names,logparams=[]): if opts.eccentricity: print(" Reading composite file, assuming eccentricity-based format ") field_names=("indx","m1", "m2", "a1x", "a1y", "a1z", "a2x", "a2y", "a2z","eccentricity", "lnL", "sigmaOverL", "ntot", "neff") +if opts.hyperbolic: + print(" Reading composite file, assuming hyperbolic-based format ") + field_names=("indx","m1", "m2", "a1x", "a1y", "a1z", "a2x", "a2y", "a2z","E0", "p_phi0", "lnL", "sigmaOverL", "ntot", "neff") field_formats = [np.float32 for x in field_names] composite_dtype = [ (x,float) for x in field_names] #np.dtype(names=field_names ,formats=field_formats) # Load posterior files diff --git a/MonteCarloMarginalizeCode/Code/bin/util_CharacterizePosterior.py b/MonteCarloMarginalizeCode/Code/bin/util_CharacterizePosterior.py index 78ae72b1..f73681e8 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_CharacterizePosterior.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_CharacterizePosterior.py @@ -36,6 +36,13 @@ dat_1d = (dat['m2']*dat['a2z']+dat['m1']*dat['a1z'])/(dat['m1']+dat['m2']) elif param=='mtot': dat_1d = dat['m1']+dat['m2'] + elif param=='b_hyp': + h = dat['E0'] + j = dat['p_phi0'] + nu = dat['m1']*dat['m2']/((dat['m1']+dat['m2'])**2) + E_eff = 1. + (h**2 - 1)/(2*nu) + b = j*h / (np.sqrt(E_eff**2 - 1)) + dat_1d = b else: dat_1d = dat[param] quant_here = np.percentile(dat_1d,100*quantile_list) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_CleanILE.py b/MonteCarloMarginalizeCode/Code/bin/util_CleanILE.py index 1377f4c2..731997e8 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_CleanILE.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_CleanILE.py @@ -56,7 +56,7 @@ indx, m1,m2, s1x,s1y,s1z,s2x,s2y,s2z,ecc, lnL, sigmaOverL, ntot, neff = line col_intrinsic = 10 if opts.hyperbolic: - indx, m1,m2, s1x,s1y,s1z,s2x,s2y,s2z,E0, p_phi0, lnL, sigmaOverL, ntot, neff = line + indx, m1, m2, s1x, s1y, s1z, s2x, s2y, s2z, E0, p_phi0, lnL, sigmaOverL, ntot, neff = line col_intrinsic = 11 elif len(line) == 13 and (not tides_on) and (not distance_on): # strip lines with the wrong length indx, m1,m2, s1x,s1y,s1z,s2x,s2y,s2z,lnL, sigmaOverL, ntot, neff = line diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py index 3ebc742e..7d230335 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ConstructIntrinsicPosterior_GenericCoordinates.py @@ -343,6 +343,9 @@ def extract_combination_from_LI(samples_LI, p): parser.add_argument("--assume-eos-but-primary-bh",action='store_true',help="Special case of known EOS, but primary is a BH") parser.add_argument("--use-eccentricity", action="store_true") parser.add_argument("--use-hyperbolic", action="store_true") +parser.add_argument("--force-scatter",default=False,action='store_true', help='For hyperbolic analyses forces only scatter grid points') +parser.add_argument("--force-plunge",default=False,action='store_true', help='For hyperbolic analyses forces only plunge grid points') +parser.add_argument("--force-zoomwhirl",default=False,action='store_true', help='For hyperbolic analyses forces only zoomwhirl grid points') parser.add_argument("--tripwire-fraction",default=0.05,type=float,help="Fraction of nmax of iterations after which n_eff needs to be greater than 1+epsilon for a small number epsilon") # FIXME hacky options added by me (Liz) to try to get my capstone project to work. @@ -362,10 +365,12 @@ def extract_combination_from_LI(samples_LI, p): opts.no_adapt_parameter =[] # needs to default to empty list ECC_MAX = opts.ecc_max ECC_MIN = opts.ecc_min + E0_MAX = opts.E0_max E0_MIN = opts.E0_min PPHI0_MAX = opts.pphi0_max -PPHI0_Min = opts.pphi0_min +PPHI0_MIN = opts.pphi0_min + no_plots = no_plots | opts.no_plots lnL_shift = 0 lnL_default_large_negative = -500 @@ -826,6 +831,12 @@ def eccentricity_prior(x): def precession_prior(x): return 0.5*np.ones(x.shape) # uniform over the interval [0.0, 2.0] +def initial_energy_prior(x): + return np.ones(x.shape) / (E0_MAX-E0_MIN) # uniform over the interval [E0_MIN, E0_MAX] + +def initial_angmom_prior(x): + return np.ones(x.shape) / (PPHI0_MAX-PPHI0_MIN) # uniform over the interval [PPHI0_MIN, PPHI0_MAX] + def unnormalized_uniform_prior(x): return np.ones(x.shape) def unnormalized_log_prior(x): @@ -873,6 +884,8 @@ def normalized_zbar_prior(z): # Other priors 'eccentricity':eccentricity_prior, 'chi_pavg':precession_prior, + 'E0':initial_energy_prior, + 'p_phi0':initial_angmom_prior, 'mu1': unnormalized_log_prior, 'mu2': unnormalized_uniform_prior } @@ -890,7 +903,9 @@ def normalized_zbar_prior(z): 'lambda_plus':[0.01,lambda_plus_max], 'lambda_minus':[-lambda_max,lambda_max], # will include the true region always...lots of overcoverage for small lambda, but adaptation will save us. 'eccentricity':[ECC_MIN, ECC_MAX], - 'chi_pavg':[0.0,2.0], + 'chi_pavg':[0.0,2.0], + 'E0':[E0_MIN,E0_MAX], + 'p_phi0':[PPHI0_MIN,PPHI0_MAX], # strongly recommend you do NOT use these as parameters! Only to insure backward compatibility with LI results 'LambdaTilde':[0.01,5000], 'DeltaLambdaTilde':[-500,500], @@ -2999,6 +3014,43 @@ def parse_corr_params(my_str): if opts.verbose: print(" Sample: Skipping " , line, ' due to ', p, val, downselect_dict[p]) + if opts.force_scatter: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-scatter points from the hyperbolic grid + hypclass = Pgrid.extract_param('hypclass') + if hypclass == 'scatter': + include_item = True + else: + include_item = False + + if opts.force_plunge: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-plunge points from the hyperbolic grid + hypclass = Pgrid.extract_param('hypclass') + if hypclass == 'plunge': + include_item = True + else: + include_item = False + + if opts.force_zoomwhirl: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-zoomwhirl points from the hyperbolic grid + hypclass = Pgrid.extract_param('hypclass') + if hypclass == 'zoomwhirl': + include_item = True + else: + include_item = False + + # Set some superfluous quantities, needed only for PN approximants, so the result is generated sensibly Pgrid.ampO =opts.amplitude_order Pgrid.phaseO =opts.phase_order diff --git a/MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py b/MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py index 371ec654..6d95abb3 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_FrameZeroNoiseSNR.py @@ -15,6 +15,7 @@ parser.add_argument("--fmin-snr", type=float, default=20) parser.add_argument("--fmax-snr", type=float, default=1500) parser.add_argument("--plot-sanity",action='store_true') +parser.add_argument("--verbose", default=False, action='store_true') opts= parser.parse_args() fminSNR = opts.fmin_snr @@ -48,7 +49,7 @@ else: analyticPSD_Q=False print("Reading PSD for instrument %s from %s" % (ifo, psd_name[ifo])) - psd_dict[ifo] = lalsimutils.load_resample_and_clean_psd(psd_name[ifo], ifo, df) + psd_dict[ifo] = lalsimutils.load_resample_and_clean_psd(psd_name[ifo], ifo, df, verbose=opts.verbose) IP = lalsimutils.ComplexIP(fLow=fminSNR, fNyq=fSample/2,deltaF=df,psd=psd_dict[ifo],fMax=fmaxSNR,analyticPSD_Q=analyticPSD_Q) rhoDet = rho_dict[ifo] = IP.norm(data_dict[ifo]) print(ifo, rho_dict[ifo]) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py b/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py index 66bae5ed..a1a0091b 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_LALWriteFrame.py @@ -40,6 +40,9 @@ parser.add_argument("--mass2",default=1.4,type=float,help='Mass 2 (solar masses)') parser.add_argument("--verbose", action="store_true",default=False) parser.add_argument("--l-max",default=4,type=float,help='Lmax number of modes') +parser.add_argument('--gen-hlmoft', action='store_true', help='Creates hoft from hlmoft') +parser.add_argument('--hyperbolic', action='store_true', help='skips tapering for hyperbolic waveforms') +parser.add_argument('--force-hyperbolic-22', action='store_true', help='Forces just the 22 modes for hyperbolic waveforms') opts= parser.parse_args() @@ -51,7 +54,8 @@ P.randomize(aligned_spin_Q=True,default_inclination=opts.incl) P.m1 = opts.mass1*lalsimutils.lsu_MSUN P.m2 = opts.mass2*lalsimutils.lsu_MSUN - P.taper = lalsimutils.lsu_TAPER_START + if not(opts.hyperbolic): + P.taper = lalsimutils.lsu_TAPER_START P.tref =1000000000 # default if opts.approx: P.approx = lalsim.GetApproximantFromString(str(opts.approx)) @@ -65,10 +69,12 @@ xmldoc = utils.load_filename(filename, verbose = True, contenthandler =lalsimutils.cthdler) sim_inspiral_table = lsctables.SimInspiralTable.get_table(xmldoc) P.copy_sim_inspiral(sim_inspiral_table[int(event)]) - P.taper = lalsimutils.lsu_TAPER_START + if not(opts.hyperbolic): + P.taper = lalsimutils.lsu_TAPER_START if opts.approx: P.approx = lalsim.GetApproximantFromString(str(opts.approx)) -P.taper = lalsimutils.lsu_TAPER_START # force taper +if not(opts.hyperbolic): + P.taper = lalsimutils.lsu_TAPER_START # force taper P.detector = opts.instrument if opts.approx == "EccentricTD": P.phaseO = 3 @@ -86,7 +92,14 @@ # Generate signal -hoft = lalsimutils.hoft(P,Lmax=opts.l_max) # include translation of source, but NOT interpolation onto regular time grid +if opts.gen_hlmoft: + if opts.hyperbolic and opts.force_hyperbolic_22: + hlmT = lalsimutils.hlmoft(P, Lmax=opts.l_max, force_22_mode=True) + else: + hlmT = lalsimutils.hlmoft(P, Lmax=opts.l_max) + hoft = lalsimutils.hoft_from_hlm(hlmT, P, return_complex=False) +else: + hoft = lalsimutils.hoft(P,Lmax=opts.l_max) # include translation of source, but NOT interpolation onto regular time grid # zero pad to be opts.seglen long, if necessary if opts.seglen/hoft.deltaT > hoft.data.length: TDlenGoal = int(opts.seglen/hoft.deltaT) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py b/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py index 42e184b6..d09ab416 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ManualOverlapGrid.py @@ -245,7 +245,15 @@ parser.add_argument("--verbose", action="store_true",default=False, help="Extra warnings") parser.add_argument("--extra-verbose", action="store_true",default=False, help="Lots of messages") parser.add_argument("--save-plots",default=False,action='store_true', help="Write plots to file (only useful for OSX, where interactive is default") +parser.add_argument('--force-scatter',default=False,action='store_true',help='For hyperbolic analyses forces only scatter grid points.') +parser.add_argument('--force-plunge',default=False,action='store_true',help='For hyperbolic analyses forces only plunge grid points.') +parser.add_argument('--force-zoomwhirl',default=False,action='store_true',help='For hyperbolic analyses forces only zoomwhirl grid points.') opts= parser.parse_args() + +force_options = [opts.force_scatter, opts.force_plunge, opts.force_zoomwhirl] # Add more if needed +if sum(bool(x) for x in force_options) > 1: + parser.error("CANNOT use multiple --force-X options at the same time!") + if opts.inj_file_out: opts.fname = opts.inj_file_out.replace(".xml.gz","") @@ -366,6 +374,43 @@ def evaluate_overlap_on_grid(hfbase,param_names, grid): for param in downselect_dict: if Pgrid.extract_param(param) < downselect_dict[param][0] or Pgrid.extract_param(param) > downselect_dict[param][1]: include_item =False + + if opts.force_scatter: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-scatter points from the hyperbolic grid + hypclass = Pgrid.extract_param('hypclass') + if hypclass == 'scatter': + include_item = True + else: + include_item = False + + if opts.force_plunge: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-plunge points from the hyperbolic grid + hypclass = Pgrid.extract_param('hypclass') + if hypclass == 'plunge': + include_item = True + else: + include_item = False + + if opts.force_zoomwhirl: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-zoomwhirl points from the hyperbolic grid + hypclass = Pgrid.extract_param('hypclass') + if hypclass == 'zoomwhirl': + include_item = True + else: + include_item = False + if include_item: grid_revised.append(line) if Pgrid.m2 <= Pgrid.m1: # do not add grid elements with m2> m1, to avoid possible code pathologies ! diff --git a/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py b/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py index 30a239b6..cb22b9a8 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_ParameterPuffball.py @@ -54,13 +54,20 @@ parser.add_argument("--downselect-parameter-range",action='append',type=str) parser.add_argument("--enforce-duration-bound",default=None,type=float,help="If present, enforce a duration bound. Used to prevent grid placement for obscenely long signals, when the window size is prescribed") parser.add_argument("--regularize",action='store_true',help="Add some ad-hoc terms based on priors, to help with nearly-singular matricies") +parser.add_argument('--force-scatter',default=False,action='store_true',help='For hyperbolic analyses forces only scatter grid points.') +parser.add_argument('--force-plunge',default=False,action='store_true',help='For hyperbolic analyses forces only plunge grid points.') +parser.add_argument('--force-zoomwhirl',default=False,action='store_true',help='For hyperbolic analyses forces only zoomwhirl grid points.') opts= parser.parse_args() +force_options = [opts.force_scatter, opts.force_plunge, opts.force_zoomwhirl] # Add more if needed +if sum(bool(x) for x in force_options) > 1: + parser.error("CANNOT use multiple --force-X options at the same time!") + if opts.random_parameter is None: opts.random_parameter = [] # Extract parameter names -coord_names = opts.parameter # Used in fit +coord_names = list(set(opts.parameter)) # Used in fit #if opts.parameter_nofit: # coord_names = coord_names + opts.parameter_nofit if coord_names is None: @@ -261,6 +268,43 @@ def test_index_distance(indx,thresh=opts.force_away): val = val/ lal.MSUN_SI if val < downselect_dict[param][0] or val > downselect_dict[param][1]: include_item =False + + if opts.force_scatter: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-scatter points from the hyperbolic grid + hypclass = P.extract_param('hypclass') + if hypclass == 'scatter': + include_item = True + else: + include_item = False + + if opts.force_plunge: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-plunge points from the hyperbolic grid + hypclass = P.extract_param('hypclass') + if hypclass == 'plunge': + include_item = True + else: + include_item = False + + if opts.force_zoomwhirl: + if include_item==False: + # no need to evaluate if the point is already downselected out + pass + else: + # removes non-zoomwhirl points from the hyperbolic grid + hypclass = P.extract_param('hypclass') + if hypclass == 'zoomwhirl': + include_item = True + else: + include_item = False + if include_item: P_out.append(P) diff --git a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py index 37acf89c..44168db2 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_RIFT_pseudo_pipe.py @@ -92,6 +92,14 @@ def retrieve_event_from_coinc(fname_coinc): event_dict["m2"] = row.mass2 event_dict["s1z"] = row.spin1z event_dict["s2z"] = row.spin2z + try: + event_dict["E0"] = row.psi3 + except: + event_dict["E0"] = 0.0 + try: + event_dict["p_phi0"] = row.beta + except: + event_dict["p_phi0"] = 0.0 try: event_dict["eccentricity"] = row.alpha4 except: @@ -220,6 +228,7 @@ def unsafe_parse_arg_string_dict(my_argstr): parser.add_argument("--limit-mc-range",default=None,type=str,help="Pass this argumen through to the helper to set the mc range") parser.add_argument("--force-mc-range",default=None,type=str,help="Pass this argumen through to the helper to set the mc range") parser.add_argument("--force-eta-range",default=None,type=str,help="Pass this argumen through to the helper to set the eta range") +parser.add_argument("--force-mtot-range",default=None,type=str,help="Pass this argument through to the helper to set the mtot range. Overrides mc parameter with mtot parameter broadly throughout the pipeline.") parser.add_argument("--force-comp-max",default=1000,type=float,help="Provde this value to override the value of component mass in CIP provided") parser.add_argument("--force-comp-min",default=1,type=float,help="Provde this value to override the value of component mass in CIP provided") parser.add_argument("--allow-subsolar", action='store_true', help="Override limits which otherwise prevent subsolar mass PE") @@ -279,8 +288,23 @@ def unsafe_parse_arg_string_dict(my_argstr): parser.add_argument("--archive-pesummary-label",default=None,help="If provided, creates a 'pesummary' directory and fills it with this run's final output at the end of the run") parser.add_argument("--archive-pesummary-event-label",default="this_event",help="Label to use on the pesummary page itself") parser.add_argument("--internal-mitigate-fd-J-frame",default="L_frame",help="L_frame|rotate, choose method to deal with ChooseFDWaveform being in wrong frame. Default is to request L frame for inputs") +parser.add_argument("--first-iteration-jumpstart",action='store_true',help="No ILE jobs the first iteration. Assumes you already have .composite files and want to get going. Particularly helpful for subdag systems") +parser.add_argument("--use-mtot-coords",action='store_true',help="Passed to the helper to configure CIP and PUFF for mtot instead of mc.") +parser.add_argument("--force-scatter-grids",action='store_true',help="Eliminates all non-scatter intrinsic points from hyperbolic grids throughout the workflow.") +parser.add_argument("--force-plunge-grids",action='store_true',help="Eliminates all non-plunge intrinsic points from hyperbolic grids throughout the workflow.") +parser.add_argument("--force-zoomwhirl-grids",action='store_true',help="Eliminates all non-zoomwhirl intrinsic points from hyperbolic grids throughout the workflow.") +parser.add_argument("--force-hyperbolic-22", action='store_true', help='Forces just the 22 modes for hyperbolic waveforms') opts= parser.parse_args() +# Ensure --assume-hyperbolic is set when using any --force-X-grids option +# Ensure only ONE of the --force-X-grids options is set +force_grids = [opts.force_scatter_grids, opts.force_plunge_grids, opts.force_zoomwhirl_grids] +if any(force_grids) and not opts.assume_hyperbolic: + parser.error("Using --force-scatter-grids, --force-plunge-grids, or --force-zoomwhirl-grids requires --assume-hyperbolic!") + +if sum(bool(x) for x in force_grids) > 1: + parser.error("CANNOT use multiple --force-X-grids options at the same time!") + if (opts.use_ini): # Attempt to lazy-parse all command line arguments from ini file @@ -487,8 +511,7 @@ def unsafe_parse_arg_string_dict(my_argstr): if opts.assume_eccentric: is_analysis_eccentric = True if opts.assume_hyperbolic: - is_analysis_hyperbolic = True - + is_analysis_hyperbolic = True dirname_run = gwid+ "_" + opts.calibration+ "_"+ opts.approx+"_fmin" + str(fmin) +"_fmin-template"+str(fmin_template) +"_lmax"+str(opts.l_max) + "_"+opts.spin_magnitude_prior if opts.online: @@ -555,6 +578,10 @@ def unsafe_parse_arg_string_dict(my_argstr): # default value for eccentricity is 0 for 'P'! Only change this value from default if eccentricity is present, do NOT want to fill it with None in particular if not(event_dict['eccentricity'] is None): P.eccentricity = event_dict["eccentricity"] + # same for the hyperbolic params + if not(event_dict['E0'] is None): + P.E0 = event_dict['E0'] + P.p_phi0 = event_dict["p_phi0"] # Write 'target_params.xml.gz' file lalsimutils.ChooseWaveformParams_array_to_xml([P], "target_params") @@ -632,6 +659,30 @@ def unsafe_parse_arg_string_dict(my_argstr): cmd += " --assume-eccentric " if is_analysis_hyperbolic: cmd += " --assume-hyperbolic " + if not(opts.force_E0_max is None): + E0_max = opts.force_E0_max + cmd += " --E0-max {} ".format(E0_max) + if not(opts.force_E0_min is None): + E0_min = opts.force_E0_min + cmd += " --E0-min {} ".format(E0_min) + if not(opts.force_pphi0_max is None): + pphi0_max = opts.force_pphi0_max + cmd += " --pphi0-max {} ".format(pphi0_max) + if not(opts.force_pphi0_min is None): + pphi0_min = opts.force_pphi0_min + cmd += " --pphi0-min {} ".format(pphi0_min) + + if opts.force_scatter_grids: + cmd += " --force-scatter-grids " + + if opts.force_plunge_grids: + cmd += " --force-plunge-grids " + + if opts.force_zoomwhirl_grids: + cmd += " --force-zoomwhirl-grids " + + if opts.force_hyperbolic_22: + cmd += " --force-hyperbolic-22 " if opts.assume_highq: cmd+= ' --assume-highq --force-grid-stretch-mc-factor 2' # the mc range, tuned to equal-mass binaries, is probably too narrow. Workaround until fixed in helper npts_it =1000 @@ -661,6 +712,10 @@ def unsafe_parse_arg_string_dict(my_argstr): cmd += " --scale-mc-range " + str(opts.scale_mc_range).replace(' ','') if not(opts.force_eta_range is None): cmd+= " --force-eta-range " + str(opts.force_eta_range).replace(' ','') +if not(opts.force_mtot_range is None): + cmd+= " --force-mtot-range " + str(opts.force_mtot_range).replace(' ','') +if opts.use_mtot_coords: + cmd+= " --use-mtot-coords " if opts.allow_subsolar: cmd += " --allow-subsolar " if opts.force_chi_max: @@ -1057,7 +1112,12 @@ def unsafe_parse_arg_string_dict(my_argstr): line += " --ecc-min {} ".format(ecc_min) if opts.assume_hyperbolic: if not(opts.internal_use_aligned_phase_coordinates): - line = line.replace('parameter mc', 'parameter mc --parameter E0 --parameter p_phi0 --use-hyperbolic') + if not(opts.use_mtot_coords): + line = line.replace('parameter mc', 'parameter mc --parameter E0 --parameter p_phi0 --use-hyperbolic') + else: + #line = line.replace('parameter mc', 'parameter mtot --parameter E0 --parameter p_phi0 --use-hyperbolic') + line = line.replace('parameter mc', 'parameter-implied mc --parameter-nofit mtot --parameter-nofit q --parameter E0 --parameter p_phi0 --use-hyperbolic') + line = line.replace('parameter delta_mc', 'parameter-implied delta_mc') else: line = line.replace('parameter-nofit mc', 'parameter-nofit mc --parameter E0 --parameter p_phi0 --use-hyperbolic') if not(opts.force_E0_max is None): @@ -1072,6 +1132,17 @@ def unsafe_parse_arg_string_dict(my_argstr): if not(opts.force_pphi0_min is None): pphi0_min = opts.force_pphi0_min line += " --pphi0-min {} ".format(pphi0_min) + + if opts.force_scatter_grids: + line += " --force-scatter " + + if opts.force_plunge_grids: + line += " --force-plunge " + + if opts.force_zoomwhirl_grids: + line += " --force-zoomwhirl " + + if not(opts.manual_extra_cip_args is None): line += " {} ".format(opts.manual_extra_cip_args) # embed with space on each side, avoid collisions line += "\n" @@ -1133,9 +1204,28 @@ def unsafe_parse_arg_string_dict(my_argstr): # puff_params += " --parameter LambdaTilde " # should already be present puff_max_it +=5 # make sure we resolve the correlations if opts.assume_eccentric: - puff_params += " --parameter eccentricity --downselect-parameter eccentricity --downselect-parameter-range '[{},{}]' ".format(ecc_min,ecc_max) + puff_params += " --parameter eccentricity --downselect-parameter eccentricity --downselect-parameter-range '[{},{}]' ".format(ecc_min,ecc_max) if opts.assume_hyperbolic: - puff_params += " --parameter E0 --downselect-parameter E0 --downselect-parameter-range '[{},{}]' --parameter E0 --downselect-parameter p_phi0 --downselect-parameter-range '[{},{}]' ".format(E0_min,E0_max,pphi0_min,pphi0_max) + puff_params += " --parameter E0 " + if not(opts.force_E0_max is None and opts.force_E0_min is None): + E0_max = opts.force_E0_max + E0_min = opts.force_E0_min + puff_params += f" --downselect-parameter E0 --downselect-parameter-range [{E0_min},{E0_max}] " + puff_params += " --parameter p_phi0 " + if not(opts.force_pphi0_max is None and opts.force_pphi0_min is None): + pphi0_max = opts.force_pphi0_max + pphi0_min = opts.force_pphi0_min + puff_params += f" --downselect-parameter p_phi0 --downselect-parameter-range [{pphi0_min},{pphi0_max}]" + + if opts.force_scatter_grids: + puff_params += ' --force-scatter ' + + if opts.force_plunge_grids: + puff_params += ' --force-plunge ' + + if opts.force_zoomwhirl_grids: + puff_params += ' --force-zoomwhirl ' + if opts.assume_highq: puff_params = puff_params.replace(' delta_mc ', ' eta ') # use natural coordinates in the high q strategy. May want to do this always puff_max_it +=3 @@ -1254,6 +1344,8 @@ def unsafe_parse_arg_string_dict(my_argstr): cmd += " --comov-distance-reweighting --comov-distance-reweighting-exe `which make_uni_comov_skymap.py` --convert-ascii2h5-exe `which convert_output_format_ascii2h5.py` " if opts.use_gauss_early: cmd += " --cip-exe-G `which util_ConstructIntrinsicPosterior_GaussianResampling.py ` " +if opts.first_iteration_jumpstart: + cmd += " --first-iteration-jumpstart " if opts.internal_use_amr: print(" AMR prototype: Using hardcoded aligned-spin settings, assembling grid, requires coinc!") cmd += " --cip-exe `which util_AMRGrid.py ` " diff --git a/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py b/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py index 072a1feb..be2ac6d5 100755 --- a/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py +++ b/MonteCarloMarginalizeCode/Code/bin/util_SimInspiralToCoinc.py @@ -6,7 +6,7 @@ import json -def ilwd_base(): +def ilwd_base(arg=None): return 0 try: from glue.ligolw import ilwd @@ -126,6 +126,8 @@ def _empty_row(obj): sngl.snr = 20. # made up, needed for some algorithms to work sngl.alpha4 = P.eccentricity + sngl.psi3 = P.E0 + sngl.beta = P.p_phi0 # add to table sngl_table.append(sngl) diff --git a/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py new file mode 100755 index 00000000..547cccf5 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/test/single_inject/single_inject_RIFT.py @@ -0,0 +1,681 @@ +#! /usr/bin/env python +# +# GOAL: single-button injection setup from ini file +# +# + +################################################################ +#### Code Description #### +################################################################ + +# This code serves to create a full RIFT injection workflow, pulling all relevant information (including waveform parameters) from a config.ini file. +# Please see the example ini file (BBH_inject_example.ini) + + +## Basic Code Process: ## + +# The code process is as follows, starting from reading in the ini file: (example: single_inject_RIFT.py -- /home/albert.einstein/path/to/config/BBH_inject_example.ini) +# Create a MDC (mock data challenge) SimInspiral table xml from the waveform parameterrs +# uses LALWriteFrame.py to create frame files for the specified detectors, stored in the signal_frames dir (also combined_frames if there is no noise) +# if specified with --use-noise, adds noise frames and dumps the combined into the combined_frames dir +# calculates SNR with util_FrameZeroNoiseSNR.py, uses network SNR as the initial guess +# writes the coinc.xml table +# runs util_RIFT_pseudo_pipe.py, which sets up the RIFT analysis workflow. In the resultant analysis directory the user needs to only submit the dag via condor_submit_dag + + +## Options: ## + +# --use-osg will allow you to run on the IGWN pool of distributed computing resources, for which you need to have the environment variables: +# SINGULARITY_RIFT_IMAGE: path to the RIFT Singularity image (container) that will run RIFT +# SINGULARITY_BASE_EXE_DIR=/usr/local/bin/ + +# --add-extrinsic will enable calculation of the extrinsic posteriors in the final iteration + +# --force-cpu will disable the GPU code pathway, this is useful for running smaller or custom jobs + +# --use-noise will add noise to the created signal frames + +# --bypass frames will skip to the util_RIFT_pseudo_pipe.py step, assuming that all the file structure is otherwise in place. +# this is useful for recreating an analysis off of existing frames without waiting for them to generate. + +# --just-frames will exit the program after the SNR report is completed. +# when combined with --bypass frames, allows for two-stage injection creation + +# --use-hyperbolic passes further options to util_RIFT_pseudo_pipe.py, enabling analysis of hyperbolic systems + + +## Other notes: ## + +# if you're using a distance marginalization table, it will take a while to create the injection. You'll see: +# IntegrationWarning: The integral is probably divergent, or slowly convergent. +# but it is working, give it a few minutes + + + +################################################################ +#### Preamble #### +################################################################ + +import numpy as np +import argparse +import os +import sys +sys.path.append(os.getcwd()) +import shutil +import ast +import glob +import re +import configparser + +from RIFT.misc.dag_utils import mkdir +import RIFT.lalsimutils as lalsimutils +import lal +import lalsimulation as lalsim + +from gwpy.timeseries import TimeSeries +from matplotlib import pyplot as plt + +# Backward compatibility +from RIFT.misc.dag_utils import which +lalapps_path2cache = which('lal_path2cache') +if lalapps_path2cache == None: + lalapps_path2cache = which('lalapps_path2cache') + +################################################################ +#### Functions & Definitions #### +################################################################ + + +def unsafe_config_get(config,args,verbose=False): + if verbose: + print(" Retrieving ", args, end=' ') + print(" Found ",eval(config.get(*args))) + return eval( config.get(*args)) + +def add_frames(channel,input_frame, noise_frame,combined_frame): + # note - this needs to be updated to work with the current install regardless of user + exe = os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(__file__)), '../pp/add_frames.py')) + if not(os.path.dirname(exe) in os.environ['PATH'].split(os.pathsep)): + print('add_frames.py not found, adding to PATH...') + os.environ['PATH'] += os.pathsep + os.path.dirname(exe) + cmd = exe + " " + channel + " " + input_frame + " " + noise_frame + " " + combined_frame + print(cmd) + os.system(cmd) + +lsu_MSUN=lal.MSUN_SI +lsu_PC=lal.PC_SI + +################################################################ +#### Options Parse #### +################################################################ + +parser = argparse.ArgumentParser() +parser.add_argument("--use-ini", default=None, type=str, required=True, help="REQUIRED. Pass ini file for parsing. Intended to reproduce lalinference_pipe functionality. Overrides most other arguments.") +parser.add_argument("--use-osg",action='store_true',help="Attempt OSG operation. Command-line level option, not at ini level to make portable results") +parser.add_argument("--add-extrinsic",action='store_true',help="Add extrinsic posterior. Corresponds to --add-extrinsic --add-extrinsic-time-resampling --batch-extrinsic for pipeline") +parser.add_argument("--force-cpu",action='store_true',help="Forces avoidance of GPUs") +parser.add_argument("--use-noise",action='store_true',help="Combine clean signal injection with fiducial noise frames. Note that you must have pycbc installed.") +parser.add_argument("--bypass-frames",action='store_true',help="Skip making mdc and frame files, use with caution") +parser.add_argument("--just-frames",action='store_true',help="Stop after making frame files, use with caution") +parser.add_argument("--use-hyperbolic",action='store_true',help="Adds hyperbolic options, requires TEOBResumS approx. NOTE: development still in progress.") +parser.add_argument("--alternate-start",action='store_true',help="Currently just passses first-iteration-jumpstart to start from CIP") +opts = parser.parse_args() + +config = configparser.ConfigParser(allow_no_value=True) #SafeConfigParser deprecated from py3.2 + +ini_path = os.path.abspath(opts.use_ini) + +# read in config +config.read(ini_path) + +bypass_frames=opts.bypass_frames +just_frames=opts.just_frames +use_noise=opts.use_noise + +dist_marg_arg = config.get('rift-pseudo-pipe','internal-marginalize-distance') +dist_marg = True if dist_marg_arg.lower() == 'true' else False + +################################################################ +#### Initialize #### +################################################################ + +# Create, go to working directory +working_dir = config.get('init','working_directory') +print(" Working dir ", working_dir) +mkdir(working_dir) +os.chdir(working_dir) +working_dir_full = os.getcwd() + +# load IFO list from config +ifos = unsafe_config_get(config,['analysis','ifos']) + +################################################################ +#### Load Injection Parameters from config into lalsimutils #### +################################################################ + +P_list = [] +P = lalsimutils.ChooseWaveformParams() + +## Waveform Model ## +approx_str = config.get('engine','approx') +P.approx=lalsim.GetApproximantFromString(approx_str) + +## Intrinsic Parameters ## +P.m1 = float(config.get('injection-parameters','m1'))*lsu_MSUN +P.m2 = float(config.get('injection-parameters','m2'))*lsu_MSUN +P.s1x = float(config.get('injection-parameters','s1x')) +P.s1y = float(config.get('injection-parameters','s1y')) +P.s1z = float(config.get('injection-parameters','s1z')) +P.s2x = float(config.get('injection-parameters','s2x')) +P.s2y = float(config.get('injection-parameters','s2y')) +P.s2z = float(config.get('injection-parameters','s2z')) +if opts.use_hyperbolic: + P.E0 = float(config.get('injection-parameters','E0')) + P.p_phi0 = float(config.get('injection-parameters','p_phi0')) +## Extrinsic Parameters ## +P.dist = float(config.get('injection-parameters','dist'))*1e6*lsu_PC +fiducial_event_time = float(config.get('injection-parameters','event_time')) # for frame addition +P.tref = fiducial_event_time +P.theta = float(config.get('injection-parameters','dec')) +P.phi = float(config.get('injection-parameters','ra')) +P.incl = float(config.get('injection-parameters','incl')) +P.psi = float(config.get('injection-parameters','psi')) +P.phiref = float(config.get('injection-parameters','phase')) +# Other waveform settings +P.fmin = float(config.get('injection-parameters','fmin')) + +P_list.append(P) + +# create mdc.xml.gz - contains all injection params +if bypass_frames: + print("Skipping mdc creation") + pass +else: + lalsimutils.ChooseWaveformParams_array_to_xml(P_list,"mdc") + + + +################################################################ +#### Create frame files and cache file #### +################################################################ + +if bypass_frames: + print("Skipping frame creation") + indx = 0 + pass +else: + # create clean signal frames # + print(" === Writing signal files === ") + + mkdir('signal_frames') + # must be compatible with duration of the noise frames + t_start = int(fiducial_event_time)-150 + t_stop = int(fiducial_event_time)+150 + + indx = 0 # need an event number, meaningless otherwise + target_subdir = 'signal_frames/event_{}'.format(indx) # this is where the signal frames will go + mkdir(working_dir_full+"/"+target_subdir) + os.chdir(working_dir_full+"/"+target_subdir) + # here we're in the signal_frames/event_0 directory. + # Loop over instruments, write frame files and cache + for ifo in ifos: + cmd = "util_LALWriteFrame.py --inj " + working_dir_full+"/mdc.xml.gz --event {} --start {} --stop {} --instrument {} --approx {}".format(indx, t_start,t_stop,ifo, approx_str) + print(cmd) + os.system(cmd) + + # Here we plot the signal frames + gwf_files = [os.path.abspath(file) for file in os.listdir(os.getcwd()) if file.endswith('.gwf')] + + for index, path in enumerate(gwf_files): + file = os.path.basename(path) + # assign correct channel name + if file.startswith('H'): + channel = 'H1:FAKE-STRAIN' + elif file.startswith('L'): + channel = 'L1:FAKE-STRAIN' + elif file.startswith('V'): + channel = 'V1:FAKE-STRAIN' + + data = TimeSeries.read(path, channel) + strain = np.array(data.value) + sample_times = np.array(data.times) + + # plot domain - don't do this for combined frames + max_amp = np.argmax(strain) + plot_domain = [sample_times[max_amp] - 1.0, sample_times[max_amp] + 0.5] + + plt.plot(sample_times, data, label=f'{channel}') + plt.xlim(plot_domain) + plt.legend(loc='best') + plt.savefig(f'{channel}_plot.png', bbox_inches='tight') + plt.close() + + if not use_noise: + # copy the contents of signal_frames/event_0 to combined_frames/event_0 + os.chdir(working_dir_full) + mkdir('combined_frames') + os.chdir(working_dir_full) + target_subdir='combined_frames/event_{}'.format(indx) + shutil.copytree(working_dir_full+'/signal_frames/event_{}'.format(indx), target_subdir) + + else: + + print(" === Using fiducial noise frames === ") + + noise_frame_dir = config.get("make_data","fiducial_noise_frames") + + print(" === Joining synthetic signals to reference noise === ") + + seglen_actual = t_stop - t_start + os.chdir(working_dir_full) + mkdir('combined_frames') + os.chdir(working_dir_full) + target_subdir='combined_frames/event_{}'.format(indx) + mkdir(target_subdir) + os.chdir(working_dir_full+"/"+target_subdir) + for ifo in ifos: + fname_input = working_dir_full+"/signal_frames/event_{}/{}-fake_strain-{}-{}.gwf".format(indx,ifo[0],int(t_start),seglen_actual) + fname_noise= noise_frame_dir + "/" + ifo + ("/%08d.gwf" % indx) + fname_output = working_dir_full+"/{}/{}-combined-{}-{}.gwf".format(target_subdir,ifo[0],int(t_start),seglen_actual) + print(ifo+":FAKE-STRAIN") + print(fname_input) + print(fname_noise) + print(fname_output) + add_frames(ifo+":FAKE-STRAIN",fname_input, fname_noise,fname_output) + + # Here we plot the combined frames + gwf_files = [os.path.abspath(file) for file in os.listdir(os.getcwd()) if file.endswith('.gwf')] + + for index, path in enumerate(gwf_files): + file = os.path.basename(path) + # assign correct channel name + if file.startswith('H'): + channel = 'H1:FAKE-STRAIN' + elif file.startswith('L'): + channel = 'L1:FAKE-STRAIN' + elif file.startswith('V'): + channel = 'V1:FAKE-STRAIN' + + data = TimeSeries.read(path, channel) + + plt.plot(data, label=f'{channel}') + plt.legend(loc='best') + plt.savefig(f'{channel}_plot.png', bbox_inches='tight') + plt.close() + + # write frame paths to cache + os.chdir(working_dir_full) + target_subdir='combined_frames/event_{}'.format(indx) + os.chdir(working_dir_full+"/"+target_subdir) + cmd = "/bin/find . -name '*.gwf' | {} > signals.cache".format(lalapps_path2cache) + os.system(cmd) + + +if not opts.use_noise: + # Calculate SNRs for zero-noise frames + target_subdir='combined_frames/event_{}'.format(indx) + target_subdir_full = working_dir_full+"/"+target_subdir + os.chdir(target_subdir_full) + if os.path.exists(target_subdir_full + '/snr-report.txt'): + print('SNR report already exists, skipping') + pass + else: + cmd = "util_FrameZeroNoiseSNR.py --cache signals.cache --psd lalsim.SimNoisePSDaLIGOZeroDetHighPower" + os.system(cmd) + with open('snr-report.txt', 'r') as file: + snr_dict = ast.literal_eval(file.read()) + file.close() + snr_guess = snr_dict['Network'] # passed to pseudo_pipe as --force-hint-snr +else: + # need a way to estimate SNR for the noise frames, skipping for now + print('Skipping SNR guess with noise frames') + pass + +if opts.just_frames: + print('Exiting after frame creation!') + sys.exit(0) + +################################################################ +#### RIFT Workflow Generation #### +################################################################ + + +# check for distance marginalization lookup table - primarily for testing +if bypass_frames: + os.chdir(working_dir_full) + if os.path.exists(working_dir_full+'/'+'distance_marginalization_lookup.npz'): + print("Using existing distance marginalization table") + dist_marg_exists = True + dist_marg_path = working_dir_full + "/distance_marginalization_lookup.npz" + else: + if dist_marg: + print("Distance marginalization table not found, it will be made during the workflow") + dist_marg_exists = False +else: + dist_marg_exists = False + + +os.chdir(working_dir_full) + +if not(bypass_frames): + # write coinc file + cmd = "util_SimInspiralToCoinc.py --sim-xml mdc.xml.gz --event {}".format(indx) + for ifo in ifos: + cmd += " --ifo {} ".format(ifo) + os.system(cmd) + +# Designate rundir +rundir = config.get('init','rundir') + +#point to signals.cache +target_file = 'combined_frames/event_{}/signals.cache'.format(indx) + +# RIFT analysis setup +cmd = 'util_RIFT_pseudo_pipe.py --use-coinc `pwd`/coinc.xml --use-ini {} --use-rundir `pwd`/{} --fake-data-cache `pwd`/{}'.format(ini_path, rundir, target_file) +if opts.add_extrinsic: + cmd += ' --add-extrinsic --add-extrinsic-time-resampling --batch-extrinsic ' +if not opts.use_noise: + cmd += ' --force-hint-snr {} '.format(snr_guess) +if dist_marg_exists: + cmd += ' --internal-marginalize-distance-file {} '.format(dist_marg_path) +if opts.force_cpu: + if opts.use_osg: + print("Don't force CPUs on the OSG, exiting.") + sys.exit(1) + else: + print('Avoiding GPUs') + cmd += ' --ile-no-gpu ' +if opts.use_osg: + print('Configuring for OSG operation') + cmd += ' --use-osg ' + cmd += ' --use-osg-file-transfer ' + cmd += ' --ile-retries 10 ' + # note: --use-osg-cip not used by default - can set this in the ini file +if opts.use_hyperbolic: + print('Adding hyperbolic options') + cmd += ' --assume-hyperbolic ' + +if opts.alternate_start: + # ensure there is first iteration jumpstart + if config.getboolean('rift-pseudo-pipe', 'first-iteration-jumpstart'): + pass + else: + cmd += ' --first-iteration-jumpstart ' + +os.system(cmd) + +################################################################ +#### Post-Workflow Cleanup #### +################################################################ + +# copy over PSDs +psd_dict = unsafe_config_get(config,["make_psd",'fiducial_psds']) +for psd in psd_dict: + psd_path = psd_dict[psd] + shutil.copy2(psd_path, rundir) + +# copy config file into rundir/reproducibility +repro_path = working_dir_full+"/"+rundir+"/"+"reproducibility" +shutil.copy2(ini_path, repro_path) + +if opts.use_osg: + # make frames_dir and copy frame files for transfer + os.chdir(working_dir_full+"/"+rundir) + mkdir('frames_dir') + os.chdir(working_dir_full+"/combined_frames/event_{}".format(indx)) + for filename in os.listdir(os.getcwd()): + if filename.endswith('.gwf'): + pathname = os.path.join(os.getcwd(), filename) + if os.path.isfile(pathname): + shutil.copy2(pathname, working_dir_full+"/"+rundir+"/frames_dir") + + +# for ease of testing, copy the distance marginalization table to the top-level dir +if dist_marg: + if os.path.exists(working_dir_full+'/'+'distance_marginalization_lookup.npz'): + pass + else: + print('Copying distance marginalization') + os.chdir(working_dir_full+"/"+rundir) + shutil.copy2('distance_marginalization_lookup.npz', working_dir_full) + + +# edit dag file + change name of dag +os.chdir(working_dir_full+"/"+rundir) +if opts.alternate_start: + with open('marginalize_intrinsic_parameters_BasicIterationWorkflow.dag', 'r') as file: + lines = file.readlines() + + # remove first 6 lines + lines = lines[6:] + + # Find the start of the parent/child section + parent_child_start = None + for i, line in enumerate(lines): + if line.startswith('PARENT'): + parent_child_start = i + break + + if parent_child_start is not None: + # Remove the first 2 lines in the parent/child section + lines = lines[:parent_child_start] + lines[parent_child_start+2:] + + # Write the modified lines back to the file + with open('marginalize_intrinsic_parameters_BasicIterationWorkflow.dag', 'w') as file: + file.writelines(lines) + + +os.rename('marginalize_intrinsic_parameters_BasicIterationWorkflow.dag', working_dir + '.dag') + +# For alternate start, copy over all.net file +if opts.alternate_start: + net_file = os.path.join(working_dir_full, 'analysis_0/all.net') + shutil.copy2(net_file, working_dir_full+"/"+rundir+"/all.net") + + +## Update ILE requirements to avoid bad hosts ## + +excluded_hosts = [] + +host_file = '/home/richard.oshaughnessy/igwn_feedback/rift_avoid_hosts.txt' + +with open(host_file, 'r') as file: + for line in file: + hostname = line.strip() + excluded_hosts.append(f"(TARGET.Machine =!= '{hostname}')") + +avoid_string = "&&".join(excluded_hosts) + +# List of file paths to modify +file_paths = ["ILE.sub","ILE_puff.sub","iteration_4_cip/ILE.sub","iteration_4_cip/ILE_puff.sub"] + +if opts.use_hyperbolic: + if config.getboolean('rift-pseudo-pipe', 'internal-use-aligned-phase-coordinates'): + file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"] + else: + file_paths = ["ILE.sub","ILE_puff.sub","iteration_3_cip/ILE.sub","iteration_3_cip/ILE_puff.sub"] + +if opts.add_extrinsic: + file_paths.append("ILE_extr.sub") + +for file_path in file_paths: + # Construct the full file path + full_path = os.path.join(os.getcwd(), file_path) + + # Read the contents of the file + with open(full_path, "r") as file: + lines = file.readlines() + + # Find the index of the line containing the 'requirements' command + requirements_index = None + for i, line in enumerate(lines): + if line.strip().startswith("requirements"): + requirements_index = i + break + + # Find the index of the line containing the 'queue' command + queue_index = len(lines) - 1 + for i, line in enumerate(reversed(lines)): + if line.strip().startswith("queue"): + queue_index = len(lines) - 1 - i + break + + # Find the index of the line containing 'request_memory' + memory_index = None + for i, line in enumerate(lines): + if line.strip().startswith("request_memory"): + memory_index = i + break + + # Modify lines + if requirements_index is not None: + #requirements_line = lines[requirements_index] + if opts.use_osg: + requirements_line = "requirements = (HAS_SINGULARITY=?=TRUE)" + "&&" + avoid_string + "\n" + lines[requirements_index] = requirements_line + # Insert the 'require_gpus' line before the 'queue' command + lines.insert(queue_index, "require_gpus = Capability >= 3.5\n") + else: + requirements_line = "requirements = " + avoid_string + "\n" + lines[requirements_index] = requirements_line + + # Modify the 'request_memory' line if opts.use_hyperbolic is true + if opts.use_hyperbolic and memory_index is not None: + lines[memory_index] = "request_memory = 16384M\n" + + # Write the modified contents back to the file + with open(full_path, "w") as file: + file.writelines(lines) + +# +# +## special configuration for hyperbolic analysis in the IGWN pool +if opts.use_hyperbolic and opts.use_osg: + + print('Modifying submit files for hyperbolics in the IGWN pool') + + # need to edit all ILE and CIP submit files + if config.getboolean('rift-pseudo-pipe', 'internal-use-aligned-phase-coordinates'): + ILE_file_paths = ["ILE.sub","ILE_puff.sub","iteration_2_cip/ILE.sub","iteration_2_cip/ILE_puff.sub"] + CIP_file_paths = ["CIP.sub", "CIP_worker.sub", "CIP_0.sub", "CIP_worker0.sub", "CIP_1.sub", "CIP_worker1.sub", "CIP_2.sub", "CIP_worker2.sub", "iteration_2_cip/CIP.sub", "iteration_2_cip/CIP_worker.sub"] + PUFF_file_paths = ["PUFF.sub", "iteration_2_cip/PUFF.sub"] + else: + ILE_file_paths = ["ILE.sub","ILE_puff.sub","iteration_3_cip/ILE.sub","iteration_3_cip/ILE_puff.sub"] + CIP_file_paths = ["CIP.sub", "CIP_worker.sub", "CIP_0.sub", "CIP_worker0.sub", "CIP_1.sub", "CIP_worker1.sub", "CIP_2.sub", "CIP_worker2.sub", "iteration_3_cip/CIP.sub", "iteration_3_cip/CIP_worker.sub"] + PUFF_file_paths = ["PUFF.sub", "iteration_3_cip/PUFF.sub"] + + + if opts.add_extrinsic: + ILE_file_paths.append("ILE_extr.sub") + + + total_file_paths = ILE_file_paths + CIP_file_paths + + + def modify_submit_file(input_file): + + #current_singularity_image = 'rift_o4b_jl-chadhenshaw-teobresums_eccentric-2024-05-14_12-02-56.sif' + #sif_path = 'osdf:///igwn/cit/staging/james.clark/' + current_singularity_image + + #current_singularity_image = 'rift_o4b_jl-2024-07-16_12-21-57.sif' + # sif_path = 'osdf:///igwn/cit/staging/chad.henshaw/' + current_singularity_image + + current_singularity_image = '/cvmfs/singularity.opensciencegrid.org/james-clark/research-projects-rit/containers-rift_o4b_jl-chadhenshaw-teobresums_eccentric:latest' + + with open(input_file, 'r') as file: + lines = file.readlines() + + #lines_to_remove = ['getenv', '+SingularityBindCVMFS', '+SingularityImage'] + lines_to_remove = ['getenv'] + #lines_to_add = [f'MY.SingularityImage = "{current_singularity_image}"', 'use_oauth_services = scitokens'] + lines_to_add = [] + #modifications = {'request_disk': 'request_disk = 3G', 'requirements': 'requirements = (HAS_SINGULARITY=?=TRUE)&&(IS_GLIDEIN=?=TRUE)', 'transfer_input_files': None} + modifications = {'request_disk': 'request_disk = 50M', 'requirements': 'requirements = (HAS_SINGULARITY=?=TRUE)&&(IS_GLIDEIN=?=TRUE)&&(TARGET.Machine =!= "deepclean.ldas.cit")'} + + modified_lines = [] + queue_line = None + + for line in lines: + stripped_line = line.strip() + + # remove lines + if any(stripped_line.startswith(keyword) for keyword in lines_to_remove): + continue + + # modify lines + modified = False + for key, new_value in modifications.items(): + if stripped_line.startswith(key): + if key == 'transfer_input_files': + #line = stripped_line + ',' + sif_path + '\n' + print('skipping this for cvfms') + else: + line = new_value + '\n' + modified_lines.append(line) + modified = True + break + + if stripped_line.startswith('queue'): + queue_line = line + elif not modified: + modified_lines.append(line) + + + # add new lines + modified_lines.extend([line + '\n' for line in lines_to_add]) + + if queue_line: + modified_lines.append(queue_line) + + with open(input_file, 'w') as file: + file.writelines(modified_lines) + + return + + def modify_CIP_coords(input_file, mtot_range): + with open(input_file, 'r') as file: + lines = file.readlines() + + # Format the new mtot_range string + mtot_range_str = f"[{mtot_range[0]}, {mtot_range[1]}]" + + for i, line in enumerate(lines): + if line.startswith("arguments ="): + # Replace --parameter mc with --parameter mtot + line = line.replace('--parameter mc', '--parameter mtot') + # Use regular expression to find and replace the --mc-range argument + line = re.sub(r'--mc-range\s*\'\[.*?\]\'', f"--mtot-range '{mtot_range_str}'", line) + lines[i] = line + + with open(input_file, 'w') as file: + file.writelines(lines) + + return + + def modify_puff_coords(input_file): + with open(input_file, 'r') as file: + lines = file.readlines() + + for i, line in enumerate(lines): + if line.startswith("arguments ="): + # Replace --parameter mc with --parameter mtot + line = line.replace('--parameter mc', '--parameter mtot') + + lines[i] = line + with open(input_file, 'w') as file: + file.writelines(lines) + + return + + for submit_file in total_file_paths: + modify_submit_file(submit_file) + + if opts.alternate_start: + # calculate mtotal range + mtot_inj = float(config.get('injection-parameters','m1')) + float(config.get('injection-parameters','m2')) + mtot_range = [0.5*mtot_inj, 1.5*mtot_inj] + for submit_file in CIP_file_paths: + modify_CIP_coords(submit_file, mtot_range) + + for submit_file in PUFF_file_paths: + modify_puff_coords(submit_file)