From 2ca0c5597590e2043bd280dd8aee3d9d09bae29a Mon Sep 17 00:00:00 2001 From: shivesh Date: Sun, 22 Apr 2018 23:18:44 -0500 Subject: Sun Apr 22 23:18:44 CDT 2018 --- bout/plot.py | 129 ++++++++++++++++++++++++++++++++++++ bout/plot_corr.py | 193 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ bout/plot_full.py | 127 +++++++++++++++++++++++++++++++++++ 3 files changed, 449 insertions(+) create mode 100644 bout/plot.py create mode 100644 bout/plot_corr.py create mode 100644 bout/plot_full.py (limited to 'bout') diff --git a/bout/plot.py b/bout/plot.py new file mode 100644 index 0000000..717cf81 --- /dev/null +++ b/bout/plot.py @@ -0,0 +1,129 @@ +import os + +import numpy as np +import numpy.ma as ma + +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import pyplot as plt +from matplotlib.offsetbox import AnchoredText +from matplotlib import rc + +rc('text', usetex=False) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +fix_sfr_mfr = [ + (1, 1, 1, 1, 2, 0), + # (1, 1, 1, 1, 0, 0), + (1, 1, 1, 0, 1, 0), +] + +# FR +# dimension = [3, 6] +dimension = [3, 6] +sigma_ratio = ['0.01'] +energy_dependance = 'spectral' +spectral_index = -2 +binning = [1e4, 1e7, 5] +fix_mixing = 'False' +fix_mixing_almost = 'False' +scale_region = "1E10" + +# Likelihood +likelihood = 'golemfit' +confidence = 2.71 # 90% for 1DOF +outformat = ['png'] + + +def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01): + mr = np.array(measured_ratio) / float(np.sum(measured_ratio)) + sr = np.array(source_ratio) / float(np.sum(source_ratio)) + si = sigma_ratio + out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format( + int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), + int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension + ) + return out + + +def get_units(dimension): + if dimension == 3: return r' / GeV' + if dimension == 4: return r'' + if dimension == 5: return r' / GeV^{-1}' + if dimension == 6: return r' / GeV^{-2}' + if dimension == 7: return r' / GeV^{-3}' + if dimension == 8: return r' / GeV^{-4}' + + +def myround(x, base=5, up=False, down=False): + if up == down and up is True: assert 0 + if up: return int(base * np.round(float(x)/base-0.5)) + elif down: return int(base * np.round(float(x)/base+0.5)) + else: int(base * np.round(float(x)/base)) + + +colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + +for i_dim, dim in enumerate(dimension): + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + yranges = [np.inf, -np.inf] + legend_handles = [] + xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$'] + ax.set_xlim(0, len(xticks)+1) + ax.set_xticklabels([''] + xticks + ['']) + ax.set_xlabel(r'BSM operator angle') + ylabel = r'${\rm log}_{10} \Lambda' + get_units(dim) + r'$' + ax.set_ylabel(ylabel) + for i_frs, frs in enumerate(fix_sfr_mfr): + print '== DIM{0}'.format(dim) + print '== FRS = {0}'.format(frs) + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index) + infile = outchain_head + '/angles_limit/fr_anfr_evidence'+ gen_identifier(frs[:3], frs[-3:], dim) + '.npy' + try: + array = np.load(infile) + except IOError: + print 'failed to open {0}'.format(infile) + continue + print 'array', array + print 'array', array.shape + for i_th in xrange(len(xticks)): + scale, llhs = array[i_th].T + min_llh = np.min(llhs) + delta_llh = 2*(llhs - min_llh) + print 'scale', scale + print 'delta_llh', delta_llh + al = scale[delta_llh < confidence] + if len(al) > 0: + label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5]) + lim = al[0] + print 'frs, dim, lim = ', frs, dim, lim + if lim < yranges[0]: yranges[0] = lim + if lim > yranges[1]: yranges[1] = lim+4 + line = plt.Line2D( + (i_th+1-0.1, i_th+1+0.1), (lim, lim), lw=3, color=colour[i_frs], label=label + ) + ax.add_line(line) + if i_th == 0: legend_handles.append(line) + x_offset = i_frs*0.05 - 0.05 + ax.annotate( + s='', xy=(i_th+1+x_offset, lim), xytext=(i_th+1+x_offset, lim+3), + arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[i_frs]} + ) + else: + print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, min_llh) + try: + yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) + # ax.set_ylim(yranges) + ax.set_ylim([-30, -20]) + except: pass + + ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right', + title='dimension {0}'.format(dim)) + for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) + for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + + for of in outformat: + fig.savefig('../images/freq/lim_DIM{0}.'.format(dim)+of, bbox_inches='tight', dpi=150) diff --git a/bout/plot_corr.py b/bout/plot_corr.py new file mode 100644 index 0000000..a75fe8a --- /dev/null +++ b/bout/plot_corr.py @@ -0,0 +1,193 @@ +import os + +import numpy as np +import numpy.ma as ma + +import scipy.interpolate as interpolate + +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import pyplot as plt +from matplotlib.offsetbox import AnchoredText +from matplotlib import rc + +rc('text', usetex=False) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +fix_sfr_mfr = [ + (1, 1, 1, 1, 2, 0), + # (1, 1, 1, 1, 0, 0), + (1, 1, 1, 0, 1, 0), +] + +# FR +# dimension = [3, 6] +dimension = [3, 6] +sigma_ratio = ['0.01'] +energy_dependance = 'spectral' +spectral_index = -2 +binning = [1e4, 1e7, 5] +fix_mixing = 'False' +fix_mixing_almost = 'False' +scale_region = "1E10" + +# Likelihood +likelihood = 'golemfit' +confidence = 4.61 # 90% for 2DOF +outformat = ['png'] + + +def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01): + mr = np.array(measured_ratio) / float(np.sum(measured_ratio)) + sr = np.array(source_ratio) / float(np.sum(source_ratio)) + si = sigma_ratio + out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format( + int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), + int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension + ) + return out + + +def get_units(dimension): + if dimension == 3: return r' / GeV' + if dimension == 4: return r'' + if dimension == 5: return r' / GeV^{-1}' + if dimension == 6: return r' / GeV^{-2}' + if dimension == 7: return r' / GeV^{-3}' + if dimension == 8: return r' / GeV^{-4}' + + +def myround(x, base=5, up=False, down=False): + if up == down and up is True: assert 0 + if up: return int(base * np.round(float(x)/base-0.5)) + elif down: return int(base * np.round(float(x)/base+0.5)) + else: int(base * np.round(float(x)/base)) + + +colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + +labels = [r'$sin^2(2\mathcal{O}_{12})$', r'$sin^2(2\mathcal{O}_{13})$', r'$sin^2(2\mathcal{O}_{23})$'] +for i_dim, dim in enumerate(dimension): + for i_frs, frs in enumerate(fix_sfr_mfr): + print '== DIM{0}'.format(dim) + print '== FRS = {0}'.format(frs) + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index) + infile = outchain_head + '/angles_corr/fr_co_evidence'+ gen_identifier(frs[:3], frs[-3:], dim) + '.npy' + # infile = '../mnrun/fr_co_evidence_033_033_033_0010_sfr_033_066_000_DIM6_single_scale.npy' + try: + array = ma.masked_invalid(np.load(infile)) + except IOError: + print 'failed to open {0}'.format(infile) + continue + print 'array', array + print 'array', array.shape + for i_scen in xrange(len(labels)): + d = array[i_scen].reshape(len(array[i_scen])**2, 3) + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + xranges = [np.inf, -np.inf] + legend_handles = [] + ax.set_ylim(0, 1) + ax.set_ylabel(labels[i_scen]) + xlabel = r'${\rm log}_{10} \Lambda' + get_units(dim) + r'$' + ax.set_xlabel(xlabel) + + x = d[:,0] + y = d[:,1] + z = d[:,2] + + print 'x', x + print 'y', y + print 'z', z + null_idx = np.argmin(z) + null = z[null_idx] + print 'null = {0}, for scale = {1}'.format(null, x[null_idx]) + z = 2*(z - null) + print 'scale', x + print 'delta_llh', z + + # x_ = np.linspace(np.min(x), np.max(x), 30) + # y_ = np.linspace(np.min(y), np.max(y), 30) + # z_ = interpolate.gridddata((x, y), z, (x_[None,:], y_[:,None]), method='nearest') + + data = np.array([x, y, z, np.ones(x.shape)]).T + print 'data', data + data_clean = [] + for d in data: + if not np.any(np.isnan(d)): data_clean.append(d) + data = np.vstack(data_clean) + sort_column = 3 + data_sorted = data[data[:,sort_column].argsort()] + uni, c = np.unique(data[:,sort_column], return_counts=True) + print uni, c + print len(uni) + print np.unique(c) + + n = len(uni) + assert len(np.unique(c)) == 1 + c = c[0] + col_array = [] + for col in data_sorted.T: + col_array.append(col.reshape(n, c)) + col_array = np.stack(col_array) + sep_arrays = [] + for x_i in xrange(n): + sep_arrays.append(col_array[:,x_i]) + + print len(sep_arrays) + sep_arrays = sep_arrays[0][:3] + print sep_arrays + + llh_90_percent = (sep_arrays[2] < confidence) + data_90_percent = sep_arrays.T[llh_90_percent].T + print 'data_90_percent', data_90_percent + + ax.tick_params(axis='x', labelsize=11) + ax.tick_params(axis='y', labelsize=11) + + mini, maxi = np.min(x), np.max(x) + ax.set_xlim((mini, maxi)) + ax.set_ylim(0, 1) + ax.grid(b=False) + + x_v = data_90_percent[0].round(decimals=4) + y_v = data_90_percent[1].round(decimals=4) + uniques = np.unique(x_v) + print 'uniques', uniques + if len(uniques) == 1: continue + bw = np.min(np.diff(uniques)) + print 'bw', bw + print np.diff(uniques) + uni_x_split = np.split(uniques, np.where(np.diff(uniques) > bw*(1e20))[0] + 1) + print 'len(uni_x_split)', len(uni_x_split) + for uni_x in uni_x_split: + p_x_l, p_y_l = [], [] + p_x_u, p_y_u = [], [] + for uni in uni_x: + idxes = np.where(x_v == uni)[0] + ymin, ymax = 1, 0 + for idx in idxes: + if y_v[idx] < ymin: ymin = y_v[idx] + if y_v[idx] > ymax: ymax = y_v[idx] + p_x_l.append(uni) + p_y_l.append(ymin) + p_x_u.append(uni) + p_y_u.append(ymax) + p_x_l, p_y_l = np.array(p_x_l, dtype=np.float64), np.array(p_y_l, dtype=np.float64) + p_x_u, p_y_u = np.array(list(reversed(p_x_u)), dtype=np.float64), np.array(list(reversed(p_y_u)), dtype=np.float64) + p_x = np.hstack([p_x_l, p_x_u]) + p_y = np.hstack([p_y_l, p_y_u]) + p_x = np.r_[p_x, p_x[0]] + p_y = np.r_[p_y, p_y[0]] + print 'p_x', p_x + print 'p_y', p_y + try: + tck, u = interpolate.splprep([p_x, p_y], s=1e-5, per=True) + xi, yi = interpolate.splev(np.linspace(0, 1, 1000), tck) + except: + xi, yi = p_x, p_y + ax.fill(xi, yi, 'r', edgecolor='r', linewidth=1) + + for of in outformat: + plt.savefig('../images/freq/lim_corr_DIM{0}_AN{1}_FRS{2}'.format(dim, i_scen, i_frs)+of, bbox_inches='tight', dpi=150) + diff --git a/bout/plot_full.py b/bout/plot_full.py new file mode 100644 index 0000000..f2e1919 --- /dev/null +++ b/bout/plot_full.py @@ -0,0 +1,127 @@ +import os + +import numpy as np +import numpy.ma as ma + +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import pyplot as plt +from matplotlib.offsetbox import AnchoredText +from matplotlib import rc + +rc('text', usetex=False) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +fix_sfr_mfr = [ + (1, 1, 1, 1, 2, 0), + # (1, 1, 1, 1, 0, 0), + (1, 1, 1, 0, 1, 0), +] + +# FR +# dimension = [3, 6] +dimension = [3, 4, 5, 6, 7, 8] +sigma_ratio = ['0.01'] +energy_dependance = 'spectral' +spectral_index = -2 +binning = [1e4, 1e7, 5] +fix_mixing = 'False' +fix_mixing_almost = 'False' +scale_region = "1E10" + +# Likelihood +likelihood = 'golemfit' +confidence = 2.71 # 90% for 1DOF +outformat = ['png'] + + +def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01): + mr = np.array(measured_ratio) / float(np.sum(measured_ratio)) + sr = np.array(source_ratio) / float(np.sum(source_ratio)) + si = sigma_ratio + out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format( + int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), + int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension + ) + return out + + +def get_units(dimension): + if dimension == 3: return r' / GeV' + if dimension == 4: return r'' + if dimension == 5: return r' / GeV^{-1}' + if dimension == 6: return r' / GeV^{-2}' + if dimension == 7: return r' / GeV^{-3}' + if dimension == 8: return r' / GeV^{-4}' + + +def myround(x, base=5, up=False, down=False): + if up == down and up is True: assert 0 + if up: return int(base * np.round(float(x)/base-0.5)) + elif down: return int(base * np.round(float(x)/base+0.5)) + else: int(base * np.round(float(x)/base)) + + +colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + +fig = plt.figure(figsize=(7, 5)) +ax = fig.add_subplot(111) + +colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} +yranges = [np.inf, -np.inf] +legend_handles = [] +ax.set_xlim(dimension[0]-1, dimension[-1]+1) +xticks = [''] + range(dimension[0], dimension[-1]+1) + [''] +ax.set_xticklabels(xticks) +ax.set_xlabel(r'BSM operator dimension ' + r'$d$') +ax.set_ylabel(r'${\rm log}_{10} \Lambda / GeV^{-d+4}$') +for i_dim, dim in enumerate(dimension): + for i_frs, frs in enumerate(fix_sfr_mfr): + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index) + infile = outchain_head + '/bayes_factor/fr_fr_evidence' + gen_identifier(frs[:3], frs[-3:], dim) + '.npy' + try: + array = np.load(infile) + except IOError: + print 'failed to open {0}'.format(infile) + continue + print 'array', array + print 'array', array.shape + scale, llhs = array.T + print 'scale min', scale[np.argmin(llhs)] + null = llhs[np.argmin(llhs)] + # null = llhs[0] + # TODO(shivesh): negative or not? + reduced_ev = 2*(llhs - null) + print 'reduced_ev', reduced_ev + al = scale[reduced_ev < confidence] + if len(al) > 0: + label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5]) + lim = al[0] + print 'frs, dim, lim = ', frs, dim, lim + if lim < yranges[0]: yranges[0] = lim + if lim > yranges[1]: yranges[1] = lim+4 + line = plt.Line2D( + (dim-0.1, dim+0.1), (lim, lim), lw=3, color=colour[i_frs], label=label + ) + ax.add_line(line) + if i_dim == 0: legend_handles.append(line) + x_offset = i_frs*0.05 - 0.05 + ax.annotate( + s='', xy=(dim+x_offset, lim), xytext=(dim+x_offset, lim+3), + arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[i_frs]} + ) + + else: + print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, null) + # print 'scales, reduced_ev', np.dstack([scale.data, reduced_ev.data]) +yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) +ax.set_ylim(yranges) + +ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right') +for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) +for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + +for of in outformat: + fig.savefig('../images/freq/full_corr.'+of, bbox_inches='tight', dpi=150) -- cgit v1.2.3