Merge branch 'temp-RIT-Tides' of...

Merge branch 'temp-RIT-Tides' of https://github.com/oshaughn/research-projects-RIT into temp-RIT-Tides

MonteCarloMarginalizeCode/Code/factored_likelihood.py

@@ -76,7 +76,7 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict,
         inv_spec_trunc_Q=False, T_spec=0., verbose=True,
          NR_group=None,NR_param=None,
         ignore_threshold=1e-4,   # dangerous for peak lnL of 25^2/2~300 : biases
-       use_external_EOB=False,nr_lookup=False,nr_lookup_valid_groups=None,no_memory=True,perturbative_extraction=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):
+       use_external_EOB=False,nr_lookup=False,nr_lookup_valid_groups=None,no_memory=True,restricted_mode_list=None,perturbative_extraction=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):
     """
     Compute < h_lm(t) | d > and < h_lm | h_l'm' >
 
@@ -298,6 +298,17 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict,
             import sys
             sys.exit(0)
 
+    # Reduce list of modes
+    if not(restricted_mode_list is None):
+            print " ====== REDUCING MODE LIST ===== "
+            mode_list_orig = hlms.keys()
+            for mode in mode_list_orig:
+                    if not mode in restricted_mode_list:
+                            del(hlms[mode])
+                            del(hlms_conj[mode])
+            print " Restricting mode list : ", restricted_mode_list, hlms.keys(), hlms_conj.keys()
+
+    # Impose the cutoff
     if not(ignore_threshold is None) and (not ROM_use_basis):
             crossTermsFiducial = ComputeModeCrossTermIP(hlms,hlms, psd_dict[detectors[0]], 
                                                         P.fmin, fMax,
@@ -305,6 +316,8 @@ def PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict,
             theWorthwhileModes =  IdentifyEffectiveModesForDetector(crossTermsFiducial, ignore_threshold, detectors)
             # Make sure worthwhile modes satisfy reflection symmetry! Do not truncate egregiously!
             theWorthwhileModes  = theWorthwhileModes.union(  set([(p,-q) for (p,q) in theWorthwhileModes]))
+            if not(restricted_mode_list is None):
+                    theWorthwhileModes = theWorthwhileModes.intersection( set(restricted_mode_list))
             print "  Worthwhile modes : ", theWorthwhileModes
             hlmsNew = {}
             hlmsConjNew = {}

MonteCarloMarginalizeCode/Code/integrate_likelihood_extrinsic

@@ -77,6 +77,7 @@ optp.add_option("--rom-integrate-intrinsic",default=False,action='store_true',he
 optp.add_option("--nr-perturbative-extraction",default=False,action='store_true')
 optp.add_option("--nr-use-provided-strain",default=False,action='store_true')
 optp.add_option("--no-memory",default=False,action='store_true', help="At present, turns off m=0 modes. Use with EXTREME caution only if requested by model developer")
+optp.add_option("--restricted-mode-list-file",default=None,help="A list of ALL modes to use in likelihood. Incredibly dangerous. Only use when comparing with models which provide restricted mode sets, or otherwise to isolate the effect of subsets of modes on the whole")
 optp.add_option("--use-external-EOB",default=False,action='store_true')
 optp.add_option("--maximize-only",default=False, action='store_true',help="After integrating, attempts to find the single best fitting point")
 optp.add_option("--dump-lnL-time-series",default=False, action='store_true',help="(requires --sim-xml) Dump lnL(t) at the injected parameters")
@@ -158,6 +159,12 @@ opts, args = optp.parse_args()
 manual_avoid_overflow_logarithm=opts.manual_logarithm_offset
 
 
+# Load in restricted mode set, if available
+restricted_mode_list=None
+if not(opts.restricted_mode_list_file is None):
+    modes =numpy.loadtxt(opts.restricted_mode_list_file,dtype=int) # columns are l m.  Must contain all. Only integers obviously
+    restricted_mode_list = [ (l,m) for l,m in modes]
+    print " RESTRICTED MODE LIST target :", restricted_mode_list
 
 intrinsic_param_names = opts.parameter
 valid_intrinsic_param_names = ['q']
@@ -465,7 +472,7 @@ rholms_intp, cross_terms, cross_terms_V,  rholms, rest=factored_likelihood.Preco
         fiducial_epoch, t_window, P, data_dict, psd_dict, opts.l_max, fmax,
         False, inv_spec_trunc_Q, T_spec,
         NR_group=NR_template_group,NR_param=NR_template_param,
-        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,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=not opts.rom_integrate_intrinsic,ROM_limit_basis_size=opts.rom_limit_basis_size_to,no_memory=opts.no_memory)
+        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,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=not opts.rom_integrate_intrinsic,ROM_limit_basis_size=opts.rom_limit_basis_size_to,no_memory=opts.no_memory,restricted_mode_list=restricted_mode_list)
 
 ### DEBUGGING CHECK: what keys are used (for rom)
 if opts.rom_use_basis:

MonteCarloMarginalizeCode/Code/util_CleanILE.py

@@ -22,6 +22,7 @@ data_at_intrinsic = {}
 my_digits=5  # safety for high-SNR BNS
 
 tides_on = False
+distance_on = False  
 col_intrinsic = 9
 
 
@@ -31,8 +32,12 @@ for fname in sys.argv[1:]:
     for line in data:
         line = np.around(line, decimals=my_digits)
         lambda1=lambda2=0
-        if len(line) == 13 and not tides_on:  # strip lines with the wrong length
+        if 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
+        elif  len(line) == 14:
+            distance_on=True
+            col_intrinsic=10
+            indx, m1,m2, s1x,s1y,s1z,s2x,s2y,s2z,dist, lnL, sigmaOverL, ntot, neff = line
         elif len(line)==15:
             tides_on  = True
             col_intrinsic =11
@@ -57,5 +62,7 @@ for key in data_at_intrinsic:
 
     if tides_on:
         print -1,  key[0],key[1], key[2], key[3],key[4], key[5],key[6], key[7], key[8],key[9], lnLmeanMinusLmax+lnLmax, sigmaNetOverL, np.sum(ntot), -1
+    elif distance_on:
+        print -1,  key[0],key[1], key[2], key[3],key[4], key[5],key[6], key[7], key[8], lnLmeanMinusLmax+lnLmax, sigmaNetOverL, np.sum(ntot), -1
     else:
         print -1,  key[0],key[1], key[2], key[3],key[4], key[5],key[6], key[7], lnLmeanMinusLmax+lnLmax, sigmaNetOverL, np.sum(ntot), -1