aboutsummaryrefslogtreecommitdiffstats
path: root/bout
diff options
context:
space:
mode:
Diffstat (limited to 'bout')
-rw-r--r--bout/plot.py129
-rw-r--r--bout/plot_corr.py193
-rw-r--r--bout/plot_full.py127
3 files changed, 449 insertions, 0 deletions
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)
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-17 17:40:52 +0000