diff options
Diffstat (limited to 'plot_bf.py')
| -rwxr-xr-x[-rw-r--r--] | plot_bf.py | 204 |
1 files changed, 183 insertions, 21 deletions
diff --git a/plot_bf.py b/plot_bf.py index faefcde..13664b9 100644..100755 --- a/plot_bf.py +++ b/plot_bf.py @@ -24,7 +24,7 @@ from matplotlib.offsetbox import AnchoredText import fr from utils import misc as misc_utils from utils.fr import normalise_fr -from utils.plot import myround, get_units +from utils.plot import bayes_factor_plot, myround, get_units rc('text', usetex=False) @@ -33,11 +33,12 @@ 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), + (1, 1, 1, 0, 1, 0), ] # FR -dimension = [3, 4, 5, 6, 7, 8] +dimension = [3, 6] +# dimension = [3, 4, 5, 6, 7, 8] sigma_ratio = ['0.01'] energy_dependance = 'spectral' spectral_index = -2 @@ -68,15 +69,16 @@ bayes_eval_bin = True # set to 'all' to run normally # Plot plot_bayes = False -plot_angles_limit = False -plot_angles_corr = True +plot_angles_limit = True +plot_angles_corr = False outformat = ['png'] -significance = np.log(10**(1/2.)) -# significance = np.log(10**(3/2.)) +# 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) @@ -96,7 +98,7 @@ for i_dim, dim in enumerate(dimension): 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'.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) + 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) @@ -110,18 +112,33 @@ for i_dim, dim in enumerate(dimension): 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 - for sc in scan_scales: + 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,:] - if plot_angles_corr: - # TODO(shivesh) - assert 0 - if len(lf.shape) == 3: lf = lf[:,0,:] raw.append(lf) except IOError: fail += 1 @@ -141,13 +158,28 @@ for i_dim, dim in enumerate(dimension): 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: - # TODO(shivesh) 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)) @@ -160,7 +192,7 @@ if plot_bayes: 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}$') + 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 @@ -187,8 +219,11 @@ if plot_bayes: else: print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, null) - yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) - ax.set_ylim(yranges) + # 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(): @@ -197,7 +232,7 @@ if plot_bayes: ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) for of in outformat: - fig.savefig('./images/bayes_factor.'+of, bbox_inches='tight', dpi=150) + 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'} @@ -210,15 +245,16 @@ if plot_angles_limit: 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'$' + 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[0] + 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] @@ -249,4 +285,130 @@ if plot_angles_limit: ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1) for of in outformat: - fig.savefig('./images/angles_limit_DIM{0}'.format(dim)+'.'+of, bbox_inches='tight', dpi=150) + 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) + |