rift_O4b_JL - hyperbolic merge request

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

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

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",

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

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: