Merge branch 'shio_banksim_dev' into 'main'

Make bank sim plots robust

See merge request !20

bin/manifold_cbc_bank_dag

@@ -55,6 +55,7 @@ with open(arg[0]) as f:
 mkdir(config['group']['out-dir'])
 mkdir(config['patches']['out-dir'])
 mkdir(config['test']['out-dir'])
+mkdir(config['test-plot']['out-dir'])
 
 #
 # Setup the seed bank to parallelize bank creation
@@ -90,7 +91,7 @@ for n in range(config['test']['num-jobs']):
     test_layer += Node(arguments = Option("key", "test"), inputs = [Option("yaml", arg), Option("bank", config['constrain']['output-h5'], remap = True)], outputs = Option("output-h5", fn, remap = True)) 
 
 test_add_layer += Node(arguments = Option("key", "test-add"), inputs = [Argument("tests", test_files, remap = True), Option("yaml", arg, remap = True)], outputs = Option("output-h5", config['test-add']['output-h5'], remap = True))
-test_plot_layer += Node(arguments = Option("key", "test-plot"), inputs = [Argument("test", config['test-add']['output-h5'], remap = True), Option("yaml", arg, remap = True)], outputs = Option("output-png", config['test-plot']['output-png'], remap = True))
+test_plot_layer += Node(arguments = Option("key", "test-plot"), inputs = [Argument("test", config['test-add']['output-h5'], remap = True), Option("yaml", arg, remap = True)])
 
 dag.attach(bank_layer)
 dag.attach(add_layer)

bin/manifold_cbc_bank_test_plot

@@ -7,7 +7,6 @@ from manifold.bin.cbc import override, test_plot
 
 # NOTE for condor file transfer to work, you need to specify input and output files here. If you do, they will override the yaml versions
 p = OptionParser(usage="Usage: %prog [options]")
-p.add_option("", "--output-png", help="Specify the output png file, will override the YAML config")
 p.add_option("", "--test", help="Specify the test h5 file, will override the YAML config")
 p.add_option("", "--key", help="Specify an optional key to start reading the config from in the YAML file")
 p.add_option("", "--yaml", help="Specify the YAML file for configuration")
@@ -16,5 +15,5 @@ p.add_option("", "--yaml", help="Specify the YAML file for configuration")
 with open(opt.yaml) as f:
     config = yaml.load(f)
     if opt.key is not None: config = config[opt.key]
-    override(config, opt, "output-png", "test")
+    override(config, opt, "test")
     test_plot(config)

examples/O4-bank-dag.yaml

@@ -69,4 +69,4 @@ test-add:
 
 test-plot:
   test: *testfinalname
-  output-png: test.png
+  out-dir: plots 

manifold/bin/cbc.py

@@ -18,6 +18,10 @@ matplotlib.rcParams.update({
     "font.family": "serif"
 })
 from matplotlib import pyplot
+from matplotlib.figure import Figure
+from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
+
+from math import ceil, floor, log10
 
 import manifold.utilities.common
 import manifold.utilities.data
@@ -25,6 +29,21 @@ from manifold import cover
 from manifold import io
 from manifold.sources import cbc
 
+def compute_mchirp(m1, m2):
+	return (m1 * m1 * m1 * m2 * m2 * m2 / (m1 + m2))**0.2
+
+def compute_eta(m1, m2):
+	return (m1 * m2)/(m1 + m2)**2.
+
+def compute_M(m1, m2):
+	return (m1 + m2)
+
+def compute_chi(m1, m2, s1z, s2z):
+	return (m1 * s1z + m2 * s2z) / (m1 + m2)
+
+def compute_mratio(m1, m2):
+	return (m1 / m2)
+
 def override(config, opt, *opts):
     for o in opts:
         thisopt = getattr(opt, o.replace("-","_"))
@@ -209,54 +228,149 @@ test: test.h5
 output-png: test.png
     """
     test = io.h5read(config['test'])
-    
+    matches = None
+
+    # obtain parameters
     mismatches = test['mismatches']
     mismatches[mismatches > 0.05] = 0.05
-    ix = numpy.argsort(mismatches)
-    mismatches =  mismatches[ix]
+    matches = 1 - mismatches
+    ix = numpy.argsort(matches)
+    matches = matches[ix]
     m1s = test['m1s'][ix]
     m2s = test['m2s'][ix]
     s1s = test['s1s'][ix]
     s2s = test['s2s'][ix]
-    chis = (m1s * s1s + m2s * s2s) / (m1s + m2s)
-    fig = pyplot.figure(figsize=(7,2.5))
-    pyplot.subplot(131)
-    pyplot.scatter(m1s, m2s, c=mismatches, s=2)
-    ax = pyplot.gca()
-    pyplot.xlabel('m1')
-    pyplot.ylabel('m2')
-    ax.set_yscale('log')
-    ax.set_xscale('log')
-    pyplot.colorbar()
-    pyplot.subplot(132)
-    pyplot.scatter(m1s, chis, c=mismatches, s=2)
-    ax = pyplot.gca()
-    pyplot.xlabel('m1')
-    pyplot.ylabel('chi')
-    ax.set_yscale('linear')
-    ax.set_xscale('log')
-    pyplot.colorbar()
-    pyplot.subplot(133)
-    mismatches.sort()
-    percent = numpy.arange(1, len(mismatches)+1) / len(mismatches)
-    mm50 = mismatches[numpy.searchsorted(percent, 0.50)]
-    mm90 = mismatches[numpy.searchsorted(percent, 0.90)]
-    mm99 = mismatches[numpy.searchsorted(percent, 0.99)]
-    mm999 = mismatches[numpy.searchsorted(percent, 0.999)]
-    pyplot.plot(mismatches, numpy.arange(len(mismatches)) / len(mismatches))
-    pyplot.text(mm50-.004,0,'50\%',rotation=90)
-    pyplot.vlines(mm50, 0, 1)
-    pyplot.text(mm90-.004,0,'90\%',rotation=90)
-    pyplot.vlines(mm90, 0, 1)
-    pyplot.text(mm99-.004,0,'99\%',rotation=90)
-    pyplot.vlines(mm99, 0, 1)
-    pyplot.text(mm999-.004,0,'99.9\%',rotation=90)
-    pyplot.vlines(mm999, 0, 1)
-    #pyplot.hist(mismatches, 25, cumulative=True)
-    pyplot.xlabel('mismatch')
-    pyplot.ylabel('fraction')
-    fig.tight_layout()
-    pyplot.savefig(config['output-png'])
+
+    # derived quatities
+    min_match = matches[int(floor(len(matches)*0.9))]
+    smallest_match = matches.min()
+    inj_M = compute_M(m1s, m2s)
+    inj_eta = compute_eta(m1s, m2s)
+    inj_mchirp = compute_mchirp(m1s, m2s)
+    inj_chi = compute_chi(m1s, m2s, s1s, s2s)
+    inj_mratio = compute_mratio(m1s, m2s)
+
+    # plotting
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    ax.plot(m1s, matches, marker = '.', linestyle="None")
+    ax.set_xlabel("Injection $m_1$ ($M_\odot$)", fontsize=12)
+    ax.set_ylabel("Fitting Factor Between Injection and Template Bank", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_injm1.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    ax.plot(m2s, matches, marker = '.', linestyle="None")
+    ax.set_xlabel("Injection $m_2$ ($M_\odot$)", fontsize=12)
+    ax.set_ylabel("Fitting Factor Between Injection and Template Bank", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_injm2.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    ax.plot(inj_M, matches, marker = ".",linestyle="None")
+    ax.set_xlabel("Injection $M$ ($M_\odot$)", fontsize=12)
+    ax.set_ylabel("Fitting Factor Between Injection and Template Bank", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_injM.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    ax.plot(inj_mchirp, matches, marker = ".",linestyle="None")
+    ax.set_xlabel("Injection $M$ ($M_\odot$)", fontsize=12)
+    ax.set_ylabel("Fitting Factor Between Injection and Template Bank", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_injmchirp.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    ax.plot(inj_eta, matches, marker = ".", linestyle="None")
+    ax.set_xlabel("Injection $\eta$", fontsize=12)
+    ax.set_ylabel("Match Between Injection and Template Bank", fontsize=12)
+    ax.set_xlim([inj_eta.min(), 0.251])
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_injeta.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    collection = ax.scatter(m1s, m2s, c=matches, s=20, vmin=smallest_match, linewidth=0, alpha=0.5, vmax=1)
+    ax.set_xlabel("Injected $m_1$ ($M_\odot$)", fontsize=12)
+    ax.set_ylabel("Injected $m_2$ ($M_\odot$)", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    fig.colorbar(collection, ax=ax).set_label("Fitting Factor", fontsize=12)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_injm2_vs_injm1.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    collection = ax.scatter(s1s, s2s, c=matches, s=20, vmin=smallest_match, linewidth=0, alpha=0.5, vmax=1)
+    ax.set_xlabel("Injected $s_{1z}$", fontsize=12)
+    ax.set_ylabel("Injected $s_{2z}$", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    fig.colorbar(collection, ax=ax).set_label("Fitting Factor", fontsize=12)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_injs1z_vs_injs2z.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    ax.plot(inj_mratio, matches, marker = ".", linestyle="None")
+    ax.set_xlabel("Injection mass ratio", fontsize=12)
+    ax.set_ylabel("Fitting Factor Between Injection and Template Bank", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_inj_mass_ratio.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_axes((0.1, 0.1, 0.85, 0.85), xscale="log", yscale="log")
+    count, bin_edges = numpy.histogram(numpy.log(1 - matches + 1e-5), bins=40) # fudge factor to avoid log(0)
+    bin_centers = numpy.exp(-(bin_edges[:-1] * bin_edges[1:])**0.5)
+    ax.plot(bin_centers, count, linewidth=2.5)
+    ax.plot((1 - min_match, 1 - min_match), (1, 10**ceil(log10(count.max()))), "k--")
+    # ax.plot((1 - min_match, 1 - min_match), (1, numpy.power(10,5)), "k--")
+    # ax.set_xlim([numpy.power(10,-5.0), numpy.power(10,0)])
+    # ax.set_ylim([1, numpy.power(10,5)])
+    ax.set_xlabel("mismatch (1 - fitting factor)", fontsize=12)
+    ax.set_ylabel("Number", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.set_title(r'$N_\mathrm{inj}=%d$, 10th percentile=%.2f' % (len(matches), min_match,), fontsize=12)
+    ax.grid(True)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_hist.png", bbox_inches='tight')
+
+    fig = Figure()
+    ax = fig.add_subplot(111)
+    collection = ax.scatter(inj_M, inj_chi, c=matches, s=20, vmin=smallest_match, linewidth=0, alpha=0.5, vmax=1)
+    ax.set_xlabel("Injected Total Mass", fontsize=12)
+    ax.set_ylabel("Injected $%s$" % "\chi_\mathrm{eff}", fontsize=12)
+    ax.tick_params(axis='x', labelsize=12)
+    ax.tick_params(axis='y', labelsize=12)
+    ax.grid(True)
+    fig.colorbar(collection, ax=ax).set_label("Fitting Factor", fontsize=12)
+    canvas = FigureCanvas(fig)
+    fig.savefig("match_vs_injM_vs_injchi.png", bbox_inches='tight')
 
 def unique_temp_ids(ids):
 	s = set()