diff options
Diffstat (limited to 'plot_bf.py')
| -rwxr-xr-x | plot_bf.py | 414 |
1 files changed, 0 insertions, 414 deletions
diff --git a/plot_bf.py b/plot_bf.py deleted file mode 100755 index 13664b9..0000000 --- a/plot_bf.py +++ /dev/null @@ -1,414 +0,0 @@ -#! /usr/bin/env python -# author : S. Mandalia -# s.p.mandalia@qmul.ac.uk -# -# date : April 14, 2018 - -""" -HESE BSM flavour ratio sensivity plotting script -""" - -from __future__ import absolute_import, division - -import os - -import numpy as np -import numpy.ma as ma - -import matplotlib as mpl -mpl.use('Agg') -from matplotlib import rc -from matplotlib import pyplot as plt -from matplotlib.offsetbox import AnchoredText - -import fr -from utils import misc as misc_utils -from utils.fr import normalise_fr -from utils.plot import bayes_factor_plot, myround, get_units - - -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' - -# Nuisance -convNorm = 1. -promptNorm = 0. -muonNorm = 1. -astroNorm = 6.9 -astroDeltaGamma = 2.5 - -# GolemFit -ast = 'p2_0' -data = 'real' - -# Bayes Factor -bayes_bins = 100 -bayes_live_points = 1000 -bayes_tolerance = 0.01 -bayes_eval_bin = True # set to 'all' to run normally - -# Plot -plot_bayes = False -plot_angles_limit = True -plot_angles_corr = False -outformat = ['png'] -# significance = np.log(10**(1/2.)) -significance = np.log(10**(3/2.)) - - -bayes_array = ma.masked_equal(np.zeros((len(dimension), len(fix_sfr_mfr), bayes_bins, 2)), 0) -angles_lim_array = np.zeros((len(dimension), len(fix_sfr_mfr), 3, bayes_bins, 2)) -angles_corr_array = np.zeros((len(dimension), len(fix_sfr_mfr), 3, bayes_bins, bayes_bins, 3)) -for i_dim, dim in enumerate(dimension): - if energy_dependance == 'mono': - outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/{2:.0E}'.format(likelihood, dim, en) - elif energy_dependance == 'spectral': - outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(likelihood, dim, spectral_index) - - bayes_output = 'None' - angles_lim_output = 'None' - angles_corr_output = 'None' - for sig in sigma_ratio: - for i_frs, frs in enumerate(fix_sfr_mfr): - outchains = outchain_head + '/fix_ifr/{0}/'.format(str(sig).replace('.', '_')) - if plot_bayes: - bayes_output = outchains + '/bayes_factor/' - if plot_angles_limit: - angles_lim_output = outchains + '/angles_limit/' - if plot_angles_corr: - angles_corr_output = outchains + '/angles_corr/' - - argstring = '--measured-ratio {0} {1} {2} --fix-source-ratio True --source-ratio {3} {4} {5} --dimension {6} --seed 24 --outfile {7} --run-mcmc False --likelihood {8} --plot-angles False --bayes-output {9} --angles-lim-output {10} --bayes-bins {11} --angles-corr-output {12}'.format(frs[0], frs[1], frs[2], frs[3], frs[4], frs[5], dim, outchains, likelihood, bayes_output, angles_lim_output, bayes_bins, angles_corr_output) - args = fr.parse_args(argstring) - fr.process_args(args) - # misc_utils.print_args(args) - - if plot_bayes: - infile = args.bayes_output+'/fr_evidence'+misc_utils.gen_identifier(args) - if plot_angles_limit: - infile = args.angles_lim_output+'/fr_an_evidence'+misc_utils.gen_identifier(args) - if plot_angles_corr: - infile = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args) - scan_scales = np.linspace( - np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.bayes_bins - ) - print 'scan_scales', scan_scales - raw = [] - fail = 0 - if plot_angles_corr: - scan_angles = np.linspace(0, 1, args.bayes_bins) - for i_sc, sc in enumerate(scan_scales): - if plot_angles_corr: - for i_an, an in enumerate(scan_angles): - idx = i_sc*args.bayes_bins + i_an - infile_s = infile + '_idx_{0}'.format(idx) - try: - lf = np.load(infile_s+'.npy') - print 'lf.shape', lf.shape - except IOError: - fail += 1 - print 'failed to open {0}'.format(infile_s) - lf = np.full((3, 1, 1, 3), np.nan) - pass - for x in xrange(len(lf)): - angles_corr_array[i_dim][i_frs][x][i_sc][i_an] = np.array(lf[x]) - continue - try: - infile_s = infile + '_scale_{0:.0E}'.format(np.power(10, sc)) - lf = np.load(infile_s+'.npy') - print lf.shape - if plot_angles_limit: - if len(lf.shape) == 3: lf = lf[:,0,:] - raw.append(lf) - except IOError: - fail += 1 - print 'failed to open {0}'.format(infile_s) - if plot_bayes: - raw.append([0, 0]) - if plot_angles_limit: - raw.append(np.zeros((3, 2))) - pass - print 'failed to open {0} files'.format(fail) - - if plot_bayes: - raw = np.vstack(raw) - if plot_angles_limit: - raw = np.vstack(raw).reshape(args.bayes_bins, 3, 2) - a = ma.masked_equal(np.zeros((3, args.bayes_bins, 2)), 0) - for i_x, x in enumerate(raw): - for i_y, y in enumerate(x): - a[i_y][i_x] = ma.masked_equal(y, 0) - if plot_angles_corr: - a = angles_corr_array[i_dim][i_frs] - a = ma.masked_invalid(a, 0) - # for i_sc in xrange(len(scan_scales)): - # for i_a in xrange(len(scan_angles)): - # try: - # bayes_factor_plot( - # a[i_sc,:,i_a,:][:,(0,2)], './mnrun/corr/test_corr_DIM{0}_FR{1}_AN{2}_SC{3}'.format(dim, i_frs, i_a, i_sc), ['png'], args - # ) - # except: pass - - if plot_bayes: - bayes_array[i_dim][i_frs] = ma.masked_equal(raw, 0) - bayes_factor_plot( bayes_array[i_dim][i_frs], './mnrun/test_full_DIM{0}_FR{1}'.format(dim, i_frs), ['png'], args - ) - - if plot_angles_limit: - angles_lim_array[i_dim][i_frs] = ma.masked_equal(a, 0) - for i_a, angle in enumerate(a): - bayes_factor_plot( - angle, './mnrun/test_angles_DIM{0}_FR{1}_AN{2}'.format(i_dim, i_frs, i_a), ['png'], args - ) - -if plot_bayes: - 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^{-1} / GeV^{-d+4}$') - for i_dim, dim in enumerate(dimension): - for i_frs, frs in enumerate(fix_sfr_mfr): - scale, evidences = bayes_array[i_dim][i_frs].T - null = evidences[np.argmin(scale)] - # TODO(shivesh): negative or not? - reduced_ev = -(evidences - null) - al = scale[reduced_ev > significance] - 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]) - try: - yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) - ax.set_ylim(yranges) - except: pass - - 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/bayes/bayes_factor.'+of, bbox_inches='tight', dpi=150) - -if plot_angles_limit: - 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^{-1}' + get_units(dim) + r'$' - ax.set_ylabel(ylabel) - for i_th in xrange(len(xticks)): - for i_frs, frs in enumerate(fix_sfr_mfr): - scale, evidences = angles_lim_array[i_dim][i_frs][i_th].T - null = evidences[np.argmin(scale)] - # TODO(shivesh): negative or not? - reduced_ev = -(evidences - null) - al = scale[reduced_ev > significance] - # print 'scales, reduced_ev', np.dstack([scale, reduced_ev]) - 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]} - ) - try: - yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) - ax.set_ylim(yranges) - 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/bayes/angles_limit_DIM{0}'.format(dim)+'.'+of, bbox_inches='tight', dpi=150) - -if plot_angles_corr: - 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) - array = angles_corr_array[i_dim][i_frs] - print 'array', array - print 'array.shape', array.shape - for i_scen in xrange(len(labels)): - d = array[i_scen].reshape(len(array[i_scen])**2, 3) - print 'd.shape', d.shape - 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^{-1}' + get_units(dim) + r'$' - ax.set_xlabel(xlabel) - - x = d[:,0] - y = d[:,1] - z = d[:,2] - - data_clean = [] - for id in d: - if not np.any(np.isnan(id)): data_clean.append(id) - d_c = np.vstack(data_clean) - - x = d_c[:,0] - y = d_c[:,1] - z = d_c[:,2] - - # print 'x', x - # print 'y', y - # print 'z', z - null_idx = np.argmin(x) - null = z[null_idx] - print 'null = {0}, for scale = {1}'.format(null, x[null_idx]) - z = -(z - null) - print 'scale', x - print 'bayes_factor', 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 - 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 - - allowed_bf = (sep_arrays[2] < significance) # Shade the excluded region - data_allowed_bf = sep_arrays.T[allowed_bf].T - print 'data_allowed_bf', data_allowed_bf - - ax.tick_params(axis='x', labelsize=11) - ax.tick_params(axis='y', labelsize=11) - - mini, maxi = np.min(scan_scales), np.max(scan_scales) - ax.set_xlim((mini, maxi)) - ax.set_ylim(0, 1) - ax.grid(b=False) - - x_v = data_allowed_bf[0].round(decimals=4) - y_v = data_allowed_bf[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/bayes/lim_corr_DIM{0}_AN{1}_FRS{2}'.format(dim, i_scen, i_frs)+of, bbox_inches='tight', dpi=150) - |