reogranise into a python package

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)