reogranise into a python package

1 files changed, 1030 insertions, 0 deletions

diff --git a/golemflavor/plot.py b/golemflavor/plot.py
new file mode 100644
index 0000000..d19a52e
--- /dev/null
+++ b/golemflavor/plot.py
@@ -0,0 +1,1030 @@
+# 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)