diff options
Diffstat (limited to 'utils/plot.py')
| -rw-r--r-- | utils/plot.py | 1030 |
1 files changed, 0 insertions, 1030 deletions
diff --git a/utils/plot.py b/utils/plot.py deleted file mode 100644 index d19a52e..0000000 --- a/utils/plot.py +++ /dev/null @@ -1,1030 +0,0 @@ -# author : S. Mandalia -# s.p.mandalia@qmul.ac.uk -# -# date : March 19, 2018 - -""" -Plotting functions for the BSM flavour ratio analysis -""" - -from __future__ import absolute_import, division - -import os -import socket -from copy import deepcopy - -import numpy as np -import numpy.ma as ma -from scipy.interpolate import splev, splprep -from scipy.ndimage.filters import gaussian_filter - -import matplotlib as mpl -import matplotlib.patches as patches -import matplotlib.gridspec as gridspec -mpl.use('Agg') - -from matplotlib import rc -from matplotlib import pyplot as plt -from matplotlib.offsetbox import AnchoredText -from matplotlib.lines import Line2D -from matplotlib.patches import Patch -from matplotlib.patches import Arrow - -tRed = list(np.array([226,101,95]) / 255.) -tBlue = list(np.array([96,149,201]) / 255.) -tGreen = list(np.array([170,196,109]) / 255.) - -import getdist -from getdist import plots, mcsamples - -import ternary -from ternary.heatmapping import polygon_generator - -import shapely.geometry as geometry - -from shapely.ops import cascaded_union, polygonize -from scipy.spatial import Delaunay - -from utils.enums import DataType, str_enum -from utils.enums import Likelihood, ParamTag, StatCateg, Texture -from utils.misc import get_units, make_dir, solve_ratio, interval -from utils.fr import angles_to_u, angles_to_fr, SCALE_BOUNDARIES - - -BAYES_K = 1. # Strong degree of belief. -# BAYES_K = 3/2. # Very strong degree of belief. -# BAYES_K = 2. # Decisive degree of belief - - -LV_ATMO_90PC_LIMITS = { - 3: (2E-24, 1E-1), - 4: (2.7E-28, 3.16E-25), - 5: (1.5E-32, 1.12E-27), - 6: (9.1E-37, 2.82E-30), - 7: (3.6E-41, 1.77E-32), - 8: (1.4E-45, 1.00E-34) -} - - -PS = 8.203E-20 # GeV^{-1} -PLANCK_SCALE = { - 5: PS, - 6: PS**2, - 7: PS**3, - 8: PS**4 -} - - -if os.path.isfile('./plot_llh/paper.mplstyle'): - plt.style.use('./plot_llh/paper.mplstyle') -elif os.environ.get('GOLEMSOURCEPATH') is not None: - plt.style.use(os.environ['GOLEMSOURCEPATH']+'/GolemFit/scripts/paper/paper.mplstyle') -if 'submitter' in socket.gethostname(): - rc('text', usetex=False) - -mpl.rcParams['text.latex.preamble'] = [ - r'\usepackage{xcolor}', - r'\usepackage{amsmath}', - r'\usepackage{amssymb}'] -mpl.rcParams['text.latex.unicode'] = True - - -def gen_figtext(args): - """Generate the figure text.""" - t = r'$' - if args.data is DataType.REAL: - t += r'\textbf{IceCube\:Preliminary}' + '$\n$' - elif args.data in [DataType.ASIMOV, DataType.REALISATION]: - t += r'{\rm\bf IceCube\:Simulation}' + '$\n$' - t += r'\rm{Injected\:composition}'+r'\:=\:({0})_\oplus'.format( - solve_ratio(args.injected_ratio).replace('_', ':') - ) + '$\n$' - t += r'{\rm Source\:composition}'+r'\:=\:({0})'.format( - solve_ratio(args.source_ratio).replace('_', ':') - ) + r'_\text{S}' - t += '$\n$' + r'{\rm Dimension}'+r' = {0}$'.format(args.dimension) - return t - - -def texture_label(x, dim): - cpt = r'c' if dim % 2 == 0 else r'a' - if x == Texture.OEU: - # return r'$\mathcal{O}_{e\mu}$' - return r'$\mathring{'+cpt+r'}_{e\mu}^{('+str(int(dim))+r')}$' - elif x == Texture.OET: - # return r'$\mathcal{O}_{e\tau}$' - return r'$\mathring{'+cpt+r'}_{\tau e}^{('+str(int(dim))+r')}$' - elif x == Texture.OUT: - # return r'$\mathcal{O}_{\mu\tau}$' - return r'$\mathring{'+cpt+r'}_{\mu\tau}^{('+str(int(dim))+r')}$' - else: - raise AssertionError - - -def cmap_discretize(cmap, N): - colors_i = np.concatenate((np.linspace(0, 1., N), (0.,0.,0.,0.))) - colors_rgba = cmap(colors_i) - indices = np.linspace(0, 1., N+1) - cdict = {} - for ki,key in enumerate(('red','green','blue')): - cdict[key] = [ (indices[i], colors_rgba[i-1,ki], colors_rgba[i,ki]) for i in xrange(N+1) ] - # Return colormap object. - return mpl.colors.LinearSegmentedColormap(cmap.name + "_%d"%N, cdict, 1024) - - -def get_limit(scales, statistic, args, mask_initial=False, return_interp=False): - max_st = np.max(statistic) - print 'scales, stat', zip(scales, statistic) - if args.stat_method is StatCateg.BAYESIAN: - if (statistic[0] - max_st) > np.log(10**(BAYES_K)): - raise AssertionError('Discovered LV!') - - try: - tck, u = splprep([scales, statistic], s=0) - except: - print 'Failed to spline' - # return None - raise - sc, st = splev(np.linspace(0, 1, 1000), tck) - - if mask_initial: - scales_rm = sc[sc >= scales[1]] - statistic_rm = st[sc >= scales[1]] - else: - scales_rm = sc - statistic_rm = st - - min_idx = np.argmin(scales) - null = statistic[min_idx] - # if np.abs(statistic_rm[0] - null) > 0.8: - # print 'Warning, null incompatible with smallest scanned scale! For ' \ - # 'DIM {0} [{1}, {2}, {3}]!'.format( - # args.dimension, *args.source_ratio - # ) - # null = statistic_rm[0] - if args.stat_method is StatCateg.BAYESIAN: - reduced_ev = -(statistic_rm - null) - print '[reduced_ev > np.log(10**(BAYES_K))]', np.sum([reduced_ev > np.log(10**(BAYES_K))]) - al = scales_rm[reduced_ev > np.log(10**(BAYES_K))] - else: - assert 0 - if len(al) == 0: - print 'No points for DIM {0} [{1}, {2}, {3}]!'.format( - args.dimension, *args.source_ratio - ) - return None - re = -(statistic-null)[scales > al[0]] - if np.sum(re < np.log(10**(BAYES_K)) - 0.1) >= 2: - print 'Warning, peaked contour does not exclude large scales! For ' \ - 'DIM {0} [{1}, {2}, {3}]!'.format( - args.dimension, *args.source_ratio - ) - return None - if np.sum(re >= np.log(10**(BAYES_K)) + 0.0) < 2: - print 'Warning, only single point above threshold! For ' \ - 'DIM {0} [{1}, {2}, {3}]!'.format( - args.dimension, *args.source_ratio - ) - return None - - if return_interp: - return (scales_rm, reduced_ev) - - # Divide by 2 to convert to standard SME coefficient - lim = al[0] - np.log10(2.) - # lim = al[0] - print 'limit = {0}'.format(lim) - return lim - - -def heatmap(data, scale, vmin=None, vmax=None, style='triangular'): - for k, v in data.items(): - data[k] = np.array(v) - style = style.lower()[0] - if style not in ["t", "h", 'd']: - raise ValueError("Heatmap style must be 'triangular', 'dual-triangular', or 'hexagonal'") - - vertices_values = polygon_generator(data, scale, style) - - vertices = [] - for v, value in vertices_values: - vertices.append(v) - return vertices - - -def get_tax(ax, scale, ax_labels, rot_ax_labels=False, fontsize=23): - ax.set_aspect('equal') - - # Boundary and Gridlines - fig, tax = ternary.figure(ax=ax, scale=scale) - - # Draw Boundary and Gridlines - tax.boundary(linewidth=2.0) - tax.gridlines(color='grey', multiple=scale/5., linewidth=0.5, alpha=0.7, ls='--') - # tax.gridlines(color='grey', multiple=scale/10., linewidth=0.2, alpha=1, ls=':') - - # Set Axis labels and Title - if rot_ax_labels: roty, rotz = (-60, 60) - else: roty, rotz = (0, 0) - tax.bottom_axis_label( - ax_labels[0], fontsize=fontsize+4, - position=(0.55, -0.10/2, 0.5), rotation=0 - ) - tax.right_axis_label( - ax_labels[1], fontsize=fontsize+4, - position=(2./5+0.1, 3./5+0.06, 0), rotation=roty - ) - tax.left_axis_label( - ax_labels[2], fontsize=fontsize+4, - position=(-0.15, 3./5+0.1, 2./5), rotation=rotz - ) - - # Remove default Matplotlib axis - tax.get_axes().axis('off') - tax.clear_matplotlib_ticks() - - # Set ticks - ticks = np.linspace(0, 1, 6) - tax.ticks(ticks=ticks, locations=ticks*scale, axis='lr', linewidth=1, - offset=0.03, fontsize=fontsize, tick_formats='%.1f') - tax.ticks(ticks=ticks, locations=ticks*scale, axis='b', linewidth=1, - offset=0.02, fontsize=fontsize, tick_formats='%.1f') - # tax.ticks() - - tax._redraw_labels() - - return tax - - -def project(p): - """Convert from flavour to cartesian.""" - a, b, c = p - x = a + b/2. - y = b * np.sqrt(3)/2. - return [x, y] - - -def project_toflavour(p, nbins): - """Convert from cartesian to flavour space.""" - x, y = p - b = y / (np.sqrt(3)/2.) - a = x - b/2. - return [a, b, nbins-a-b] - - -def tax_fill(ax, points, nbins, **kwargs): - pol = np.array(map(project, points)) - ax.fill(pol.T[0]*nbins, pol.T[1]*nbins, **kwargs) - - -def alpha_shape(points, alpha): - """ - Compute the alpha shape (concave hull) of a set - of points. - @param points: Iterable container of points. - @param alpha: alpha value to influence the - gooeyness of the border. Smaller numbers - don't fall inward as much as larger numbers. - Too large, and you lose everything! - """ - if len(points) < 4: - # When you have a triangle, there is no sense - # in computing an alpha shape. - return geometry.MultiPoint(list(points)).convex_hull - def add_edge(edges, edge_points, coords, i, j): - """ - Add a line between the i-th and j-th points, - if not in the list already - """ - if (i, j) in edges or (j, i) in edges: - # already added - return - edges.add( (i, j) ) - edge_points.append(coords[ [i, j] ]) - coords = np.array([point.coords[0] - for point in points]) - tri = Delaunay(coords) - edges = set() - edge_points = [] - # loop over triangles: - # ia, ib, ic = indices of corner points of the - # triangle - for ia, ib, ic in tri.vertices: - pa = coords[ia] - pb = coords[ib] - pc = coords[ic] - # Lengths of sides of triangle - a = np.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2) - b = np.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2) - c = np.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2) - # Semiperimeter of triangle - s = (a + b + c)/2.0 - # Area of triangle by Heron's formula - area = np.sqrt(s*(s-a)*(s-b)*(s-c)) - circum_r = a*b*c/(4.0*area) - # Here's the radius filter. - #print circum_r - if circum_r < 1.0/alpha: - add_edge(edges, edge_points, coords, ia, ib) - add_edge(edges, edge_points, coords, ib, ic) - add_edge(edges, edge_points, coords, ic, ia) - m = geometry.MultiLineString(edge_points) - triangles = list(polygonize(m)) - return cascaded_union(triangles), edge_points - - -def flavour_contour(frs, nbins, coverage, ax=None, smoothing=0.4, - hist_smooth=0.05, plot=True, fill=False, oversample=1., - delaunay=False, d_alpha=1.5, d_gauss=0.08, debug=False, - **kwargs): - """Plot the flavour contour for a specified coverage.""" - # Histogram in flavour space - os_nbins = nbins * oversample - H, b = np.histogramdd( - (frs[:,0], frs[:,1], frs[:,2]), - bins=(os_nbins+1, os_nbins+1, os_nbins+1), - range=((0, 1), (0, 1), (0, 1)) - ) - H = H / np.sum(H) - - # 3D smoothing - H_s = gaussian_filter(H, sigma=hist_smooth) - - # Finding coverage - H_r = np.ravel(H_s) - H_rs = np.argsort(H_r)[::-1] - H_crs = np.cumsum(H_r[H_rs]) - thres = np.searchsorted(H_crs, coverage/100.) - mask_r = np.zeros(H_r.shape) - mask_r[H_rs[:thres]] = 1 - mask = mask_r.reshape(H_s.shape) - - # Get vertices inside covered region - binx = np.linspace(0, 1, os_nbins+1) - interp_dict = {} - for i in xrange(len(binx)): - for j in xrange(len(binx)): - for k in xrange(len(binx)): - if mask[i][j][k] == 1: - interp_dict[(i, j, k)] = H_s[i, j, k] - vertices = np.array(heatmap(interp_dict, os_nbins)) - points = vertices.reshape((len(vertices)*3, 2)) - if debug: - ax.scatter(*(points/float(oversample)).T, marker='o', s=3, alpha=1.0, color=kwargs['color'], zorder=9) - - pc = geometry.MultiPoint(points) - if not delaunay: - # Convex full to find points forming exterior bound - polygon = pc.convex_hull - ex_cor = np.array(list(polygon.exterior.coords)) - else: - # Delaunay - concave_hull, edge_points = alpha_shape(pc, alpha=d_alpha) - polygon = geometry.Polygon(concave_hull.buffer(1)) - if d_gauss == 0.: - ex_cor = np.array(list(polygon.exterior.coords)) - else: - ex_cor = gaussian_filter( - np.array(list(polygon.exterior.coords)), sigma=d_gauss - ) - - # Join points with a spline - tck, u = splprep([ex_cor.T[0], ex_cor.T[1]], s=0., per=1, k=1) - xi, yi = map(np.array, splev(np.linspace(0, 1, 300), tck)) - - # Spline again to smooth - if smoothing != 0: - tck, u = splprep([xi, yi], s=smoothing, per=1, k=3) - xi, yi = map(np.array, splev(np.linspace(0, 1, 600), tck)) - - xi /= float(oversample) - yi /= float(oversample) - ev_polygon = np.dstack((xi, yi))[0] - - # Remove points interpolated outside flavour triangle - f_ev_polygon = np.array(map(lambda x: project_toflavour(x, nbins), ev_polygon)) - - xf, yf, zf = f_ev_polygon.T - mask = np.array((xf < 0) | (yf < 0) | (zf < 0) | (xf > nbins) | - (yf > nbins) | (zf > nbins)) - ev_polygon = np.dstack((xi[~mask], yi[~mask]))[0] - - # Plot - if plot: - if fill: - ax.fill( - ev_polygon.T[0], ev_polygon.T[1], - label=r'{0}\%'.format(int(coverage)), **kwargs - ) - else: - ax.plot( - ev_polygon.T[0], ev_polygon.T[1], - label=r'{0}\%'.format(int(coverage)), **kwargs - ) - else: - return ev_polygon - -def plot_Tchain(Tchain, axes_labels, ranges): - """Plot the Tchain using getdist.""" - Tsample = mcsamples.MCSamples( - samples=Tchain, labels=axes_labels, ranges=ranges - ) - - Tsample.updateSettings({'contours': [0.90, 0.99]}) - Tsample.num_bins_2D=10 - Tsample.fine_bins_2D=50 - Tsample.smooth_scale_2D=0.05 - - g = plots.getSubplotPlotter() - g.settings.num_plot_contours = 2 - g.settings.axes_fontsize = 10 - g.settings.figure_legend_frame = False - g.settings.lab_fontsize = 20 - g.triangle_plot( - [Tsample], filled=True,# contour_colors=['green', 'lightgreen'] - ) - return g - - -def chainer_plot(infile, outfile, outformat, args, llh_paramset, fig_text=None, - labels=None, ranges=None): - """Make the triangle plot.""" - if hasattr(args, 'plot_elements'): - if not args.plot_angles and not args.plot_elements: - return - elif not args.plot_angles: - return - - if not isinstance(infile, np.ndarray): - raw = np.load(infile) - else: - raw = infile - print 'raw.shape', raw.shape - print 'raw', raw - - make_dir(outfile), make_dir - if fig_text is None: - fig_text = gen_figtext(args) - - if labels is None: axes_labels = llh_paramset.labels - else: axes_labels = labels - if ranges is None: ranges = llh_paramset.ranges - - if args.plot_angles: - print "Making triangle plots" - Tchain = raw - g = plot_Tchain(Tchain, axes_labels, ranges) - - mpl.pyplot.figtext(0.6, 0.7, fig_text, fontsize=20) - - # for i_ax_1, ax_1 in enumerate(g.subplots): - # for i_ax_2, ax_2 in enumerate(ax_1): - # if i_ax_1 == i_ax_2: - # itv = interval(Tchain[:,i_ax_1], percentile=90.) - # for l in itv: - # ax_2.axvline(l, color='gray', ls='--') - # ax_2.set_title(r'${0:.2f}_{{{1:.2f}}}^{{+{2:.2f}}}$'.format( - # itv[1], itv[0]-itv[1], itv[2]-itv[1] - # ), fontsize=10) - - # if not args.fix_mixing: - # sc_index = llh_paramset.from_tag(ParamTag.SCALE, index=True) - # itv = interval(Tchain[:,sc_index], percentile=90.) - # mpl.pyplot.figtext( - # 0.5, 0.3, 'Scale 90% Interval = [1E{0}, 1E{1}], Center = ' - # '1E{2}'.format(itv[0], itv[2], itv[1]) - # ) - - for of in outformat: - print 'Saving', outfile+'_angles.'+of - g.export(outfile+'_angles.'+of) - - if not hasattr(args, 'plot_elements'): - return - - if args.plot_elements: - print "Making triangle plots" - if args.fix_mixing_almost: - raise NotImplementedError - nu_index = llh_paramset.from_tag(ParamTag.NUISANCE, index=True) - fr_index = llh_paramset.from_tag(ParamTag.MMANGLES, index=True) - sc_index = llh_paramset.from_tag(ParamTag.SCALE, index=True) - if not args.fix_source_ratio: - sr_index = llh_paramset.from_tag(ParamTag.SRCANGLES, index=True) - - nu_elements = raw[:,nu_index] - fr_elements = np.array(map(flat_angles_to_u, raw[:,fr_index])) - sc_elements = raw[:,sc_index] - if not args.fix_source_ratio: - sr_elements = np.array(map(angles_to_fr, raw[:,sr_index])) - if args.fix_source_ratio: - Tchain = np.column_stack( - [nu_elements, fr_elements, sc_elements] - ) - else: - Tchain = np.column_stack( - [nu_elements, fr_elements, sc_elements, sr_elements] - ) - - trns_ranges = np.array(ranges)[nu_index,].tolist() - trns_axes_labels = np.array(axes_labels)[nu_index,].tolist() - if args.fix_mixing is not MixingScenario.NONE: - trns_axes_labels += \ - [r'\mid \tilde{U}_{e1} \mid' , r'\mid \tilde{U}_{e2} \mid' , r'\mid \tilde{U}_{e3} \mid' , \ - r'\mid \tilde{U}_{\mu1} \mid' , r'\mid \tilde{U}_{\mu2} \mid' , r'\mid \tilde{U}_{\mu3} \mid' , \ - r'\mid \tilde{U}_{\tau1} \mid' , r'\mid \tilde{U}_{\tau2} \mid' , r'\mid \tilde{U}_{\tau3} \mid'] - trns_ranges += [(0, 1)] * 9 - if not args.fix_scale: - trns_axes_labels += [np.array(axes_labels)[sc_index].tolist()] - trns_ranges += [np.array(ranges)[sc_index].tolist()] - if not args.fix_source_ratio: - trns_axes_labels += [r'\phi_e', r'\phi_\mu', r'\phi_\tau'] - trns_ranges += [(0, 1)] * 3 - - g = plot_Tchain(Tchain, trns_axes_labels, trns_ranges) - - if args.data is DataType.REAL: - plt.text(0.8, 0.7, 'IceCube Preliminary', color='red', fontsize=15, - ha='center', va='center') - elif args.data in [DataType.ASIMOV, DataType.REALISATION]: - plt.text(0.8, 0.7, 'IceCube Simulation', color='red', fontsize=15, - ha='center', va='center') - - mpl.pyplot.figtext(0.5, 0.7, fig_text, fontsize=15) - for of in outformat: - print 'Saving', outfile+'_elements'+of - g.export(outfile+'_elements.'+of) - - -def plot_statistic(data, outfile, outformat, args, scale_param, label=None): - """Make MultiNest factor or LLH value plot.""" - print 'Making Statistic plot' - fig_text = gen_figtext(args) - if label is not None: fig_text += '\n' + label - - print 'data', data - print 'data.shape', data.shape - print 'outfile', outfile - try: - scales, statistic = ma.compress_rows(data).T - lim = get_limit(deepcopy(scales), deepcopy(statistic), args, mask_initial=True) - tck, u = splprep([scales, statistic], s=0) - except: - return - sc, st = splev(np.linspace(0, 1, 1000), tck) - scales_rm = sc[sc >= scales[1]] - statistic_rm = st[sc >= scales[1]] - - min_idx = np.argmin(scales) - null = statistic[min_idx] - # fig_text += '\nnull lnZ = {0:.2f}'.format(null) - - if args.stat_method is StatCateg.BAYESIAN: - reduced_ev = -(statistic_rm - null) - elif args.stat_method is StatCateg.FREQUENTIST: - reduced_ev = -2*(statistic_rm - null) - - fig = plt.figure(figsize=(7, 5)) - ax = fig.add_subplot(111) - - xlims = SCALE_BOUNDARIES[args.dimension] - ax.set_xlim(xlims) - ax.set_xlabel(r'${\rm log}_{10}\left[\Lambda^{-1}_{'+ \ - r'{0}'.format(args.dimension)+r'}'+ \ - get_units(args.dimension)+r'\right]$', fontsize=16) - - if args.stat_method is StatCateg.BAYESIAN: - ax.set_ylabel(r'$\text{Bayes\:Factor}\:\left[\text{ln}\left(B_{0/1}\right)\right]$') - elif args.stat_method is StatCateg.FREQUENTIST: - ax.set_ylabel(r'$-2\Delta {\rm LLH}$') - - # ymin = np.round(np.min(reduced_ev) - 1.5) - # ymax = np.round(np.max(reduced_ev) + 1.5) - # ax.set_ylim((ymin, ymax)) - - ax.scatter(scales[1:], -(statistic[1:]-null), color='r') - ax.plot(scales_rm, reduced_ev, color='k', linewidth=1, alpha=1, ls='-') - - ax.axhline(y=np.log(10**(BAYES_K)), color='red', alpha=1., linewidth=1.2, ls='--') - ax.axvline(x=lim, color='red', alpha=1., linewidth=1.2, ls='--') - - at = AnchoredText( - fig_text, prop=dict(size=10), frameon=True, loc=4 - ) - at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5") - ax.add_artist(at) - make_dir(outfile) - for of in outformat: - print 'Saving as {0}'.format(outfile+'.'+of) - fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) - - -def plot_table_sens(data, outfile, outformat, args, show_lvatmo=True): - print 'Making TABLE sensitivity plot' - argsc = deepcopy(args) - - dims = args.dimensions - srcs = args.source_ratios - if args.texture is Texture.NONE: - textures = [Texture.OET, Texture.OUT] - else: - textures = [args.texture] - - if len(srcs) > 3: - raise NotImplementedError - - xlims = (-60, -20) - ylims = (0.5, 1.5) - - colour = {2:'red', 1:'blue', 0:'green'} - rgb_co = {2:[1,0,0], 1:[0,0,1], 0:[0.0, 0.5019607843137255, 0.0]} - - fig = plt.figure(figsize=(8, 6)) - gs = gridspec.GridSpec(len(dims), 1) - gs.update(hspace=0.15) - - first_ax = None - legend_log = [] - legend_elements = [] - - for idim, dim in enumerate(dims): - print '|||| DIM = {0}'.format(dim) - argsc.dimension = dim - gs0 = gridspec.GridSpecFromSubplotSpec( - len(textures), 1, subplot_spec=gs[idim], hspace=0 - ) - - for itex, tex in enumerate(textures): - argsc.texture = tex - ylabel = texture_label(tex, dim) - # if angles == 2 and ian == 0: continue - ax = fig.add_subplot(gs0[itex]) - - if first_ax is None: - first_ax = ax - - ax.set_xlim(xlims) - ax.set_ylim(ylims) - ax.set_yticks([], minor=True) - ax.set_yticks([1.], minor=False) - ax.set_yticklabels([ylabel], fontsize=13) - ax.yaxis.tick_right() - for xmaj in ax.xaxis.get_majorticklocs(): - ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.2, linewidth=1) - ax.get_xaxis().set_visible(False) - if itex == len(textures) - 2: - ax.spines['bottom'].set_alpha(0.6) - elif itex == len(textures) - 1: - ax.text( - -0.04, 1.0, '$d = {0}$'.format(dim), fontsize=16, - rotation='90', verticalalignment='center', - transform=ax.transAxes - ) - # dim_label_flag = False - ax.spines['top'].set_alpha(0.6) - ax.spines['bottom'].set_alpha(0.6) - - for isrc, src in enumerate(srcs): - print '== src', src - argsc.source_ratio = src - - if dim in PLANCK_SCALE.iterkeys(): - ps = np.log10(PLANCK_SCALE[dim]) - if ps < xlims[0]: - ax.annotate( - s='', xy=(xlims[0], 1), xytext=(xlims[0]+1, 1), - arrowprops={'arrowstyle': '->, head_length=0.2', - 'lw': 1, 'color':'purple'} - ) - elif ps > xlims[1]: - ax.annotate( - s='', xy=(xlims[1]-1, 1), xytext=(xlims[1], 1), - arrowprops={'arrowstyle': '<-, head_length=0.2', - 'lw': 1, 'color':'purple'} - ) - else: - ax.axvline(x=ps, color='purple', alpha=1., linewidth=1.5) - - try: - scales, statistic = ma.compress_rows(data[idim][isrc][itex]).T - except: continue - lim = get_limit(deepcopy(scales), deepcopy(statistic), argsc, mask_initial=True) - if lim is None: continue - - ax.axvline(x=lim, color=colour[isrc], alpha=1., linewidth=1.5) - ax.add_patch(patches.Rectangle( - (lim, ylims[0]), 100, np.diff(ylims), fill=True, - facecolor=colour[isrc], alpha=0.3, linewidth=0 - )) - - if isrc not in legend_log: - legend_log.append(isrc) - label = r'$\left('+r'{0}'.format(solve_ratio(src)).replace('_',':')+ \ - r'\right)_{\text{S}}\:\text{at\:source}$' - legend_elements.append( - Patch(facecolor=rgb_co[isrc]+[0.3], - edgecolor=rgb_co[isrc]+[1], label=label) - ) - - if itex == len(textures)-1 and show_lvatmo: - LV_lim = np.log10(LV_ATMO_90PC_LIMITS[dim]) - ax.add_patch(patches.Rectangle( - (LV_lim[1], ylims[0]), LV_lim[0]-LV_lim[1], np.diff(ylims), - fill=False, hatch='\\\\' - )) - - ax.get_xaxis().set_visible(True) - ax.set_xlabel(r'${\rm New\:Physics\:Scale}\:[\:{\rm log}_{10} (\Lambda^{-1}_{d}\:/\:{\rm GeV}^{-d+4})\: ]$', - labelpad=5, fontsize=19) - ax.tick_params(axis='x', labelsize=16) - - purple = [0.5019607843137255, 0.0, 0.5019607843137255] - if show_lvatmo: - legend_elements.append( - Patch(facecolor='none', hatch='\\\\', edgecolor='k', label='IceCube, Nature.Phy.14,961(2018)') - ) - legend_elements.append( - Patch(facecolor=purple+[0.7], edgecolor=purple+[1], label='Planck Scale Expectation') - ) - legend = first_ax.legend( - handles=legend_elements, prop=dict(size=11), loc='upper left', - title='Excluded regions', framealpha=1., edgecolor='black', - frameon=True - ) - first_ax.set_zorder(10) - plt.setp(legend.get_title(), fontsize='11') - legend.get_frame().set_linestyle('-') - - ybound = 0.595 - if args.data is DataType.REAL: - # fig.text(0.295, 0.684, 'IceCube Preliminary', color='red', fontsize=13, - fig.text(0.278, ybound, r'\bf IceCube Preliminary', color='red', fontsize=13, - ha='center', va='center', zorder=11) - elif args.data is DataType.REALISATION: - fig.text(0.278, ybound-0.05, r'\bf IceCube Simulation', color='red', fontsize=13, - ha='center', va='center', zorder=11) - else: - fig.text(0.278, ybound, r'\bf IceCube Simulation', color='red', fontsize=13, - ha='center', va='center', zorder=11) - - make_dir(outfile) - for of in outformat: - print 'Saving plot as {0}'.format(outfile+'.'+of) - fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) - - -def plot_x(data, outfile, outformat, args, normalise=False): - """Limit plot as a function of the source flavour ratio for each operator - texture.""" - print 'Making X sensitivity plot' - dim = args.dimension - if dim < 5: normalise = False - srcs = args.source_ratios - x_arr = np.array([i[0] for i in srcs]) - if args.texture is Texture.NONE: - textures = [Texture.OEU, Texture.OET, Texture.OUT] - else: - textures = [args.texture] - - # Rearrange data structure - r_data = np.full(( - data.shape[1], data.shape[0], data.shape[2], data.shape[3] - ), np.nan) - - for isrc in xrange(data.shape[0]): - for itex in xrange(data.shape[1]): - r_data[itex][isrc] = data[isrc][itex] - r_data = ma.masked_invalid(r_data) - print r_data.shape, 'r_data.shape' - - fig = plt.figure(figsize=(7, 6)) - ax = fig.add_subplot(111) - - ylims = SCALE_BOUNDARIES[dim] - if normalise: - if dim == 5: ylims = (-24, -8) - if dim == 6: ylims = (-12, 8) - if dim == 7: ylims = (0, 20) - if dim == 8: ylims = (12, 36) - else: - if dim == 3: ylims = (-28, -22) - if dim == 4: ylims = (-35, -26) - if dim == 5: SCALE_BOUNDARIES[5] - xlims = (0, 1) - - colour = {0:'red', 2:'blue', 1:'green'} - rgb_co = {0:[1,0,0], 2:[0,0,1], 1:[0.0, 0.5019607843137255, 0.0]} - - legend_log = [] - legend_elements = [] - labelsize = 13 - largesize = 17 - - ax.set_xlim(xlims) - ax.set_ylim(ylims) - xticks = [0, 1/3., 0.5, 2/3., 1] - # xlabels = [r'$0$', r'$\frac{1}{3}$', r'$\frac{1}{2}$', r'$\frac{2}{3}$', r'$1$'] - xlabels = [r'$0$', r'$1 / 3$', r'$1/2$', r'$2/3$', r'$1$'] - ax.set_xticks([], minor=True) - ax.set_xticks(xticks, minor=False) - ax.set_xticklabels(xlabels, fontsize=largesize) - if dim != 4 or dim != 3: - yticks = range(ylims[0], ylims[1], 2) + [ylims[1]] - ax.set_yticks(yticks, minor=False) - if dim == 3 or dim == 4: - yticks = range(ylims[0], ylims[1], 1) + [ylims[1]] - ax.set_yticks(yticks, minor=False) - # for ymaj in ax.yaxis.get_majorticklocs(): - # ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.2, linewidth=1) - for xmaj in xticks: - if xmaj == 1/3.: - ax.axvline(x=xmaj, ls='--', color='gray', alpha=0.5, linewidth=0.3) - # else: - # ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.2, linewidth=1) - - ax.text( - (1/3.)+0.01, 0.01, r'$(0.33:0.66:0)_\text{S}$', fontsize=labelsize, - transform=ax.transAxes, rotation='vertical', va='bottom' - ) - ax.text( - 0.96, 0.01, r'$(1:0:0)_\text{S}$', fontsize=labelsize, - transform=ax.transAxes, rotation='vertical', va='bottom', ha='left' - ) - ax.text( - 0.01, 0.01, r'$(0:1:0)_\text{S}$', fontsize=labelsize, - transform=ax.transAxes, rotation='vertical', va='bottom' - ) - yl = 0.55 - if dim == 3: yl = 0.65 - ax.text( - 0.03, yl, r'${\rm \bf Excluded}$', fontsize=largesize, - transform=ax.transAxes, color = 'g', rotation='vertical', zorder=10 - ) - ax.text( - 0.95, 0.55, r'${\rm \bf Excluded}$', fontsize=largesize, - transform=ax.transAxes, color = 'b', rotation='vertical', zorder=10 - ) - - for itex, tex in enumerate(textures): - print '|||| TEX = {0}'.format(tex) - lims = np.full(len(srcs), np.nan) - - for isrc, src in enumerate(srcs): - x = src[0] - print '|||| X = {0}'.format(x) - args.source_ratio = src - d = r_data[itex][isrc] - if np.sum(d.mask) > 2: continue - scales, statistic = ma.compress_rows(d).T - lim = get_limit(deepcopy(scales), deepcopy(statistic), args, mask_initial=True) - if lim is None: continue - if normalise: - lim -= np.log10(PLANCK_SCALE[dim]) - lims[isrc] = lim - - lims = ma.masked_invalid(lims) - size = np.sum(~lims.mask) - if size == 0: continue - - print 'x_arr, lims', zip(x_arr, lims) - if normalise: - zeropoint = 100 - else: - zeropoint = 0 - lims[lims.mask] = zeropoint - - l0 = np.argwhere(lims == zeropoint)[0] - h0 = len(lims) - np.argwhere(np.flip(lims, 0) == zeropoint)[0] - lims[int(l0):int(h0)] = zeropoint - - x_arr_a = [x_arr[0]-0.1] + list(x_arr) - x_arr_a = list(x_arr_a) + [x_arr_a[-1]+0.1] - lims = [lims[0]] + list(lims) - lims = list(lims) + [lims[-1]] - - s = 0.2 - g = 2 - if dim == 3 and tex == Texture.OUT: - s = 0.4 - g = 4 - if dim in (4,5) and tex == Texture.OUT: - s = 0.5 - g = 5 - if dim == 7 and tex == Texture.OET: - s = 1.6 - g = 2 - if dim == 7 and tex == Texture.OUT: - s = 2.0 - g = 20 - if dim == 8 and tex == Texture.OET: - s = 0.8 - g = 6 - if dim == 8 and tex == Texture.OUT: - s = 1.7 - g = 8 - - # ax.scatter(x_arr_a, lims, color='black', s=1) - tck, u = splprep([x_arr_a, lims], s=0, k=1) - x, y = splev(np.linspace(0, 1, 200), tck) - tck, u = splprep([x, y], s=s) - x, y = splev(np.linspace(0, 1, 400), tck) - y = gaussian_filter(y, sigma=g) - ax.fill_between(x, y, zeropoint, color=rgb_co[itex]+[0.3]) - # ax.scatter(x, y, color='black', s=1) - # ax.scatter(x_arr_a, lims, color=rgb_co[itex], s=8) - - if itex not in legend_log: - legend_log.append(itex) - # label = texture_label(tex, dim)[:-1] + r'\:{\rm\:texture}$' - label = texture_label(tex, dim)[:-1] + r'\:({\rm this\:work})$' - legend_elements.append( - Patch(facecolor=rgb_co[itex]+[0.3], - edgecolor=rgb_co[itex]+[1], label=label) - ) - - LV_lim = np.log10(LV_ATMO_90PC_LIMITS[dim]) - if normalise: - LV_lim -= np.log10(PLANCK_SCALE[dim]) - ax.add_patch(patches.Rectangle( - (xlims[0], LV_lim[1]), np.diff(xlims), LV_lim[0]-LV_lim[1], - fill=False, hatch='\\\\' - )) - - if dim in PLANCK_SCALE: - ps = np.log10(PLANCK_SCALE[dim]) - if normalise and dim == 6: - ps -= np.log10(PLANCK_SCALE[dim]) - ax.add_patch(Arrow( - 0.24, -0.009, 0, -5, width=0.12, capstyle='butt', - facecolor='purple', fill=True, alpha=0.8, - edgecolor='darkmagenta' - )) - ax.add_patch(Arrow( - 0.78, -0.009, 0, -5, width=0.12, capstyle='butt', - facecolor='purple', fill=True, alpha=0.8, - edgecolor='darkmagenta' - )) - - ax.text( - 0.26, 0.5, r'${\rm \bf Quantum\:Gravity\:Frontier}$', - fontsize=largesize-2, transform=ax.transAxes, va='top', - ha='left', color='purple' - ) - if dim > 5: - ax.axhline(y=ps, color='purple', alpha=1., linewidth=1.5) - - cpt = r'c' if dim % 2 == 0 else r'a' - if normalise: - ft = r'${\rm New\:Physics\:Scale}\:[\:{\rm log}_{10} \left (\mathring{'+cpt+r'}^{(' + \ - r'{0}'.format(args.dimension)+r')}\cdot{\rm E}_{\:\rm P}' - if dim > 5: ft += r'^{\:'+ r'{0}'.format(args.dimension-4)+ r'}' - ft += r'\right )\: ]$' - fig.text( - 0.01, 0.5, ft, ha='left', - va='center', rotation='vertical', fontsize=largesize - ) - else: - fig.text( - 0.01, 0.5, - r'${\rm New\:Physics\:Scale}\:[\:{\rm log}_{10} \left (\mathring{'+cpt+r'}^{(' + - r'{0}'.format(args.dimension)+r')}\:' + get_units(args.dimension) + - r'\right )\: ]$', ha='left', - va='center', rotation='vertical', fontsize=largesize - ) - - ax.set_xlabel( - r'${\rm Source\:Composition}\:[\:\left (\:x:1-x:0\:\right )_\text{S}\:]$', - labelpad=10, fontsize=largesize - ) - ax.tick_params(axis='x', labelsize=largesize-1) - - purple = [0.5019607843137255, 0.0, 0.5019607843137255] - # legend_elements.append( - # Patch(facecolor=purple+[0.7], edgecolor=purple+[1], label='Planck Scale Expectation') - # ) - legend_elements.append( - Patch(facecolor='none', hatch='\\\\', edgecolor='k', label='IceCube [TODO]') - ) - legend = ax.legend( - handles=legend_elements, prop=dict(size=labelsize-2), - loc='upper center', title='Excluded regions', framealpha=1., - edgecolor='black', frameon=True, bbox_to_anchor=(0.5, 1) - ) - plt.setp(legend.get_title(), fontsize=labelsize) - legend.get_frame().set_linestyle('-') - - # ybound = 0.14 - # if args.data is DataType.REAL: - # fig.text(0.7, ybound, r'\bf IceCube Preliminary', color='red', fontsize=13, - # ha='center', va='center', zorder=11) - # elif args.data is DataType.REALISATION: - # fig.text(0.7, ybound-0.05, r'\bf IceCube Simulation', color='red', fontsize=13, - # ha='center', va='center', zorder=11) - # else: - # fig.text(0.7, ybound, r'\bf IceCube Simulation', color='red', fontsize=13, - # ha='center', va='center', zorder=11) - - make_dir(outfile) - for of in outformat: - print 'Saving plot as {0}'.format(outfile + '.' + of) - fig.savefig(outfile + '.' + of, bbox_inches='tight', dpi=150) |