Applied Michael's patch #4639

Original: 3aec0441b3d320ee52bb3a56d22487a9c1003330

lalinference/python/lalinference/lalinference_pipe_utils.py

@@ -489,7 +489,7 @@ def get_roq_mchirp_priors(path, roq_paths, roq_params, key, gid=None,sim_inspira
   ## XML and GID cannot be given at the same time
   ## sim_inspiral must already point at the right row
   mc_priors = {}
- 
+
   if gid is not None and sim_inspiral is not None:
     print "Error in get_roq_mchirp_priors, cannot use both gid and sim_inspiral\n"
     sys.exit(1)
@@ -518,6 +518,32 @@ def get_roq_mchirp_priors(path, roq_paths, roq_params, key, gid=None,sim_inspira
 
   return mc_priors, trigger_mchirp
 
+def get_roq_component_mass_priors(path, roq_paths, roq_params, key, gid=None,sim_inspiral=None):
+
+  ## XML and GID cannot be given at the same time
+  ## sim_inspiral must already point at the right row
+  m1_priors = {}
+  m2_priors = {}
+
+  if gid is not None and sim_inspiral is not None:
+    print "Error in get_roq_mchirp_priors, cannot use both gid and sim_inspiral\n"
+    sys.exit(1)
+
+  for roq in roq_paths:
+    params=os.path.join(path,roq,'params.dat')
+    roq_params[roq]=np.genfromtxt(params,names=True)
+    m1_priors[roq]=[float(roq_params[roq]['mass1min']),float(roq_params[roq]['mass1max'])]
+    m2_priors[roq]=[float(roq_params[roq]['mass2min']),float(roq_params[roq]['mass2max'])]
+
+  if gid is not None:
+    trigger_mchirp = get_trigger_chirpmass(gid)
+  elif sim_inspiral is not None:
+    trigger_mchirp = sim_inspiral.mchirp
+  else:
+    trigger_mchirp = None
+
+  return m1_priors, m2_priors, trigger_mchirp
+
 def get_roq_mass_freq_scale_factor(mc_priors, trigger_mchirp, force_flow=None):
   mc_max = mc_priors['4s'][1]
   mc_min = mc_priors['128s'][0]
@@ -534,6 +560,27 @@ def get_roq_mass_freq_scale_factor(mc_priors, trigger_mchirp, force_flow=None):
 def create_pfn_tuple(filename,protocol='file://',site='local'):
     return( (os.path.basename(filename),protocol+os.path.abspath(filename),site) )
 
+def mchirp_from_components(m1, m2):
+  return (m1*m2)**(3.0/5.0) / (m1+m2)**(1.0/5.0)
+
+def Query_ROQ_Bounds_Type(path, roq_paths):
+  # Assume that parametrization of ROQ bounds is independent of seglen; just look at first one
+  import numpy as np
+  roq = roq_paths[0]
+  params = os.path.join(path,roq,'params.dat')
+  roq_params0 = np.genfromtxt(params,names=True)
+  roq_names_set = set(roq_params0.dtype.names)
+  component_mass_bounds_set = set(['mass1min', 'mass1max', 'mass2min', 'mass2max'])
+  chirp_mass_q_bounds_set = set(['chirpmassmin', 'chirpmassmax', 'qmin', 'qmax'])
+  if roq_names_set.issuperset(component_mass_bounds_set):
+    roq_bounds = 'component_mass'
+  elif roq_names_set.issuperset(chirp_mass_q_bounds_set):
+    roq_bounds = 'chirp_mass_q'
+  else:
+    print 'Invalid bounds for ROQ. Ether (m1,m2) or (mc,q) bounds are supported.'
+    sys.exit(1)
+  return roq_bounds
+
 class LALInferencePipelineDAG(pipeline.CondorDAG):
   def __init__(self,cp,dax=False,first_dag=True,previous_dag=None,site='local'):
     self.subfiles=[]

lalinference/python/lalinference_pipe.py

@@ -94,12 +94,19 @@ if cp.has_option('paths','roq_b_matrix_directory'):
   roq_paths=os.listdir(path)
   roq_params={}
   roq_force_flow = None
+  gid=None
+  row=None
   def key(item): # to order the ROQ bases
     return float(item[1]['seglen'])
 
-  print "WARNING: Overwriting user choice of flow, srate, seglen,mc_min, mc_max and q-min"
+  if cp.has_option('lalinference','roq_force_flow'):
+     roq_force_flow = cp.getfloat('lalinference','roq_force_flow')
+     print "WARNING: Forcing the f_low to ", str(roq_force_flow), "Hz"
+
+  print "WARNING: Overwriting user choice of flow, srate, seglen, and (mc_min, mc_max and q-min) or (mass1_min, mass1_max, mass2_min, mass2_max)"
+
   if opts.gid is not None:
-    mc_priors, trigger_mchirp = pipe_utils.get_roq_mchirp_priors(path, roq_paths, roq_params, key, gid=opts.gid)
+    gid=opts.gid
   elif opts.injections is not None or cp.has_option('input','injection-file'):
     print "Only 0-th event in the XML table will be considered while running with ROQ\n"
     # Read event 0 from  Siminspiral Table
@@ -109,33 +116,42 @@ if cp.has_option('paths','roq_b_matrix_directory'):
     else:
       inxml=cp.get('input','injection-file')
     injTable=SimInspiralUtils.ReadSimInspiralFromFiles([inxml])
-    row=injTable[0] 
-
-    mc_priors, trigger_mchirp = pipe_utils.get_roq_mchirp_priors(path, roq_paths, roq_params, key, sim_inspiral=row)
-
-  if cp.has_option('lalinference','roq_force_flow'):
-     roq_force_flow = cp.getfloat('lalinference','roq_force_flow')
-     print "WARNING: Forcing the f_low to", str(roq_force_flow), "Hz"
-
-  if opts.gid is not None or (opts.injections is not None or cp.has_option('input','injection-file')): 
-        roq_mass_freq_scale_factor = pipe_utils.get_roq_mass_freq_scale_factor(mc_priors, trigger_mchirp, roq_force_flow)
-
-	for mc_prior in mc_priors:
-		mc_priors[mc_prior] = array(mc_priors[mc_prior])*roq_mass_freq_scale_factor
-	# find mass bin containing the trigger
-	trigger_bin = None
-	for roq in roq_paths:
-		print 	
-		if mc_priors[roq][0] <= trigger_mchirp <= mc_priors[roq][1]:
-			trigger_bin = roq
-			break
-	roq_paths = [trigger_bin]
+    row=injTable[0]
+  else:
+    gid=None
+    row=None
+
+  roq_bounds = pipe_utils.Query_ROQ_Bounds_Type(path, roq_paths)
+  if roq_bounds == 'chirp_mass_q':
+    print 'ROQ has bounds in chirp mass and mass-ratio'
+    mc_priors, trigger_mchirp = pipe_utils.get_roq_mchirp_priors(path, roq_paths, roq_params, key, gid=gid, sim_inspiral=row)
+  elif roq_bounds == 'component_mass':
+    print 'ROQ has bounds in component masses'
+    # get component mass bounds, then compute the chirp mass that can be safely covered
+    # further below we pass along the component mass bounds to the sampler, not the tighter chirp-mass, q bounds
+    m1_priors, m2_priors, trigger_mchirp = pipe_utils.get_roq_component_mass_priors(path, roq_paths, roq_params, key, gid=gid, sim_inspiral=row)
+    mc_priors = {}
+    for (roq,m1_prior), (roq2,m2_prior) in zip(m1_priors.items(), m2_priors.items()):
+      mc_priors[roq] = sorted([pipe_utils.mchirp_from_components(m1_prior[1], m2_prior[0]), pipe_utils.mchirp_from_components(m1_prior[0], m2_prior[1])])
+
+  roq_mass_freq_scale_factor = pipe_utils.get_roq_mass_freq_scale_factor(mc_priors, trigger_mchirp, roq_force_flow)
+
+  if opts.gid is not None or (opts.injections is not None or cp.has_option('input','injection-file')):
+
+    for mc_prior in mc_priors:
+      mc_priors[mc_prior] = array(mc_priors[mc_prior])*roq_mass_freq_scale_factor
+    # find mass bin containing the trigger
+    trigger_bin = None
+    for roq in roq_paths:
+      print
+      if mc_priors[roq][0] <= trigger_mchirp <= mc_priors[roq][1]:
+        trigger_bin = roq
+        break
+    roq_paths = [trigger_bin]
   else:
-	mc_priors, trigger_mchirp = pipe_utils.get_roq_mchirp_priors(path, roq_paths, roq_params, key, None)
-       roq_mass_freq_scale_factor = pipe_utils.get_roq_mass_freq_scale_factor(mc_priors, trigger_mchirp, roq_force_flow)
-       for mc_prior in mc_priors:
-		mc_priors[mc_prior] = array(mc_priors[mc_prior])*roq_mass_freq_scale_factor
- 
+    for mc_prior in mc_priors:
+      mc_priors[mc_prior] = array(mc_priors[mc_prior])*roq_mass_freq_scale_factor
+
 else:
   roq_paths=[None]
 
@@ -188,13 +204,11 @@ for sampler in samps:
         path=cp.get('paths','roq_b_matrix_directory')
         thispath=os.path.join(path,roq)
         cp.set('paths','roq_b_matrix_directory',thispath)
-        mc_min=mc_priors[roq][0]
-        mc_max=mc_priors[roq][1]
         flow=int(roq_params[roq]['flow'] / roq_mass_freq_scale_factor)
         srate=int(2.*roq_params[roq]['fhigh'] / roq_mass_freq_scale_factor)
 	if srate > 8192:
 		srate = 8192
-	
+
         seglen=int(roq_params[roq]['seglen'] * roq_mass_freq_scale_factor)
         # params.dat uses the convention q>1 so our q_min is the inverse of their qmax
         q_min=1./float(roq_params[roq]['qmax'])
@@ -205,11 +219,25 @@ for sampler in samps:
         for i in tmp.keys():
           tmp[i]=flow
         cp.set('lalinference','flow',str(tmp))
-        cp.set('engine','chirpmass-min',str(mc_min))
-        cp.set('engine','chirpmass-max',str(mc_max))
-        cp.set('engine','q-min',str(q_min))
-        cp.set('engine','comp-min', str(np.max(roq_params[roq]['compmin'] * roq_mass_freq_scale_factor, mc_min * np.power(1+q_min, 1./5.) * np.power(q_min, 2./5.)))
-        cp.set('engine','comp-max', str(mc_max * np.power(1+q_min, 1./5.) * np.power(q_min, -3./5.)))
+        if roq_bounds == 'chirp_mass_q':
+          mc_min=mc_priors[roq][0]
+          mc_max=mc_priors[roq][1]
+          # params.dat uses the convention q>1 so our q_min is the inverse of their qmax
+          q_min=1./float(roq_params[roq]['qmax'])
+          cp.set('engine','chirpmass-min',str(mc_min))
+          cp.set('engine','chirpmass-max',str(mc_max))
+          cp.set('engine','q-min',str(q_min))
+          cp.set('engine','comp-min', str(np.max(roq_params[roq]['compmin'] * roq_mass_freq_scale_factor, mc_min * np.power(1+q_min, 1./5.) * np.power(q_min, 2./5.)))
+          cp.set('engine','comp-max', str(mc_max * np.power(1+q_min, 1./5.) * np.power(q_min, -3./5.)))
+        elif roq_bounds == 'component_mass':
+          m1_min = m1_priors[roq][0]
+          m1_max = m1_priors[roq][1]
+          m2_min = m2_priors[roq][0]
+          m2_max = m2_priors[roq][1]
+          cp.set('engine','mass1-min',str(m1_min))
+          cp.set('engine','mass1-max',str(m1_max))
+          cp.set('engine','mass2-min',str(m2_min))
+          cp.set('engine','mass2-max',str(m2_max))
 
       if opts.run_path is not None:
         cp.set('paths','basedir',os.path.abspath(opts.run_path))