Sun May 6 21:22:46 CDT 2018

7 files changed, 349 insertions, 97 deletions

diff --git a/utils/fr.py b/utils/fr.py
index 61b5a61..5acd6f8 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -212,7 +212,7 @@ def estimate_scale(args):
         upper_s = (m_eign/args.binning[0]) / (args.binning[0]**(args.dimension-3))
         scale = np.power(10, np.average(np.log10([lower_s, upper_s])))
         diff = upper_s / lower_s
-        scale_region = (lower_s/diff, upper_s*diff)
+        scale_region = (lower_s/np.power(10, args.dimension), upper_s*diff*np.power(10, args.dimension))
         scale_region = [np.power(10, np.round(np.log10(x))) for x in scale_region]
     return scale, scale_region
 
diff --git a/utils/gf.py b/utils/gf.py
index ea4f61f..098faa9 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -57,9 +57,14 @@ def steering_params(args):
     steering_categ = args.ast
     params = gf.SteeringParams()
     params.quiet = False
+    params.fastmode = True
     params.simToLoad= steering_categ.name.lower()
     params.evalThreads = args.threads
     # params.evalThreads = thread_factors(args.threads)[1]
+
+    # For Tianlu
+    # params.years = [999]
+
     return params
 
 
diff --git a/utils/likelihood.py b/utils/likelihood.py
index 2e4a22d..c2822a1 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -167,9 +167,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
     for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)):
         param.value = flavour_angles[idx]
 
-    print 'llh_paramset', llh_paramset
-    print 'hypo_paramset', hypo_paramset
-    print 'fr', fr
     if args.likelihood is Likelihood.FLAT:
         llh = 1.
     elif args.likelihood is Likelihood.GAUSSIAN:
@@ -179,7 +176,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
         llh = gf_utils.get_llh(fitter, hypo_paramset)
     elif args.likelihood is Likelihood.GF_FREQ:
         lhh = gf_utils.get_llh_freq(fitter, hypo_paramset)
-    print 'llh', llh
     return llh
 
 
diff --git a/utils/misc.py b/utils/misc.py
index cad03bc..5483675 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -12,6 +12,7 @@ from __future__ import absolute_import, division
 import os
 import errno
 import multiprocessing
+from fractions import gcd
 
 import argparse
 from operator import attrgetter
@@ -33,14 +34,18 @@ class SortingHelpFormatter(argparse.HelpFormatter):
         super(SortingHelpFormatter, self).add_arguments(actions)
 
 
+def solve_ratio(fr):
+    denominator = reduce(gcd, fr)
+    return [int(x/denominator) for x in fr]
+
+
 def gen_identifier(args):
     f = '_DIM{0}'.format(args.dimension)
-    mr1, mr2, mr3 = args.measured_ratio
+    mr1, mr2, mr3 = solve_ratio(args.measured_ratio)
     if args.fix_source_ratio:
-        sr1, sr2, sr3 = args.source_ratio
-        f += '_sfr_{0:03d}_{1:03d}_{2:03d}_mfr_{3:03d}_{4:03d}_{5:03d}'.format(
-            int(sr1*100), int(sr2*100), int(sr3*100),
-            int(mr1*100), int(mr2*100), int(mr3*100)
+        sr1, sr2, sr3 = solve_ratio(args.source_ratio)
+        f += '_sfr_{0:G}_{1:G}_{2:G}_mfr_{3:G}_{4:G}_{5:G}'.format(
+            sr1, sr2, sr3, mr1, mr2, mr3
         )
         if args.fix_mixing:
             f += '_fix_mixing'
@@ -128,6 +133,21 @@ def make_dir(outfile):
         else:
             raise
 
+def remove_option(parser, arg):
+    for action in parser._actions:
+        if (vars(action)['option_strings']
+            and vars(action)['option_strings'][0] == arg) \
+                or vars(action)['dest'] == arg:
+            parser._remove_action(action)
+
+    for action in parser._action_groups:
+        vars_action = vars(action)
+        var_group_actions = vars_action['_group_actions']
+        for x in var_group_actions:
+            if x.dest == arg:
+                var_group_actions.remove(x)
+                return
+
 
 def seed_parse(s):
     if s.lower() == 'none':
diff --git a/utils/multinest.py b/utils/multinest.py
index 9dd0742..f3595f9 100644
--- a/utils/multinest.py
+++ b/utils/multinest.py
@@ -16,7 +16,7 @@ import numpy as np
 from pymultinest import analyse, run
 
 from utils import likelihood
-from utils.misc import gen_outfile_name, make_dir
+from utils.misc import gen_identifier, make_dir
 
 
 def CubePrior(cube, ndim, n_params):
@@ -37,13 +37,15 @@ def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
         llh_paramset[pm].value = mn_paramset[pm].value
     theta = llh_paramset.values
     # print 'llh_paramset', llh_paramset
-    return likelihood.ln_prob(
+    llh = likelihood.ln_prob(
         theta=theta,
         args=args,
         asimov_paramset=asimov_paramset,
         llh_paramset=llh_paramset,
         fitter=fitter
     )
+    # print 'llh', llh
+    return llh
 
 
 def mn_argparse(parser):
@@ -77,8 +79,8 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
         fitter          = fitter
     )
 
-    prefix = './mnrun/DIM{0}/{1}_{2:>019}_'.format(
-        args.dimension, gen_outfile_name(args), np.random.randint(0, 2**63)
+    prefix = './mnrun/DIM{0}/{1}_{2:>010}_'.format(
+        args.dimension, gen_identifier(args), np.random.randint(0, 2**30)
     )
     make_dir(prefix)
     print 'Running evidence calculation for {0}'.format(prefix)
diff --git a/utils/param.py b/utils/param.py
index f861264..25558c0 100644
--- a/utils/param.py
+++ b/utils/param.py
@@ -252,10 +252,16 @@ def get_paramsets(args, nuisance_paramset):
         ])
     if not args.fix_scale:
         tag = ParamTag.SCALE
-        llh_paramset.append(
-            Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
-                  tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag)
-        )
+        if hasattr(args, 'dimension'):
+            llh_paramset.append(
+                Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
+                      tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag)
+            )
+        elif hasattr(args, 'dimensions'):
+            llh_paramset.append(
+                Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
+                      tex=r'{\rm log}_{10}\Lambda^{-1} / GeV^{-d+4}', tag=tag)
+            )
     if not args.fix_source_ratio:
         tag = ParamTag.SRCANGLES
         llh_paramset.extend([
diff --git a/utils/plot.py b/utils/plot.py
index 3d94cc1..8f4eba8 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -10,6 +10,7 @@ Plotting functions for the BSM flavour ratio analysis
 from __future__ import absolute_import, division
 
 import os
+from copy import deepcopy
 
 import numpy as np
 import matplotlib as mpl
@@ -112,18 +113,18 @@ def plot_Tchain(Tchain, axes_labels, ranges):
 def gen_figtext(args):
     """Generate the figure text."""
     t = ''
-    mr1, mr2, mr3 = args.measured_ratio
+    mr1, mr2, mr3 = misc_utils.solve_ratio(args.measured_ratio)
     if args.fix_source_ratio:
-        sr1, sr2, sr3 = args.source_ratio
-        t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \
-                'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \
-                '{5:.2f}]\nDimension = {6}'.format(
+        sr1, sr2, sr3 = misc_utils.solve_ratio(args.source_ratio)
+        t += 'Source flavour ratio = [{0}, {1}, {2}]\nIC ' \
+                'observed flavour ratio = [{3}, {4}, ' \
+                '{5}]\nDimension = {6}'.format(
                     sr1, sr2, sr3, mr1, mr2, mr3, args.dimension,
                     int(args.energy)
                 )
     else:
-        t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \
-                '{2:.2f}]\nDimension = {3}'.format(
+        t += 'IC observed flavour ratio = [{0}, {1}, ' \
+                '{2}]\nDimension = {3}'.format(
                     mr1, mr2, mr3, args.dimension, int(args.energy)
                 )
     if args.fix_scale:
@@ -229,28 +230,33 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset):
             g.export(outfile+'_elements.'+of)
 
 
+def myround(x, base=5, up=False, down=False):
+    if up == down and up is True: assert 0
+    if up: return int(base * np.round(float(x)/base-0.5))
+    elif down: return int(base * np.round(float(x)/base+0.5))
+    else: int(base * np.round(float(x)/base))
+
+
 def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
     """Make MultiNest factor or LLH value plot."""
-    print "Making Statistic 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
     scales, statistic = data.T
+    min_idx = np.argmin(scales)
+    null = statistic[min_idx]
     if args.stat_method is StatCateg.BAYESIAN:
-        min_idx = np.argmin(scales)
-        null = statistic[min_idx]
         reduced_ev = -(statistic - null)
     elif args.stat_method is StatCateg.FREQUENTIST:
-        min_idx = np.argmin(scales)
-        null = statistic[min_idx]
         reduced_ev = -2*(statistic - null)
 
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
 
-    ax.set_xlim(scale_param.ranges)
+    ax.set_xlim(np.log10(args.scale_region))
     ax.set_xlabel('$'+scale_param.tex+'$')
     if args.stat_method is StatCateg.BAYESIAN:
         ax.set_ylabel(r'Bayes Factor')
@@ -271,85 +277,302 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
     ax.add_artist(at)
     misc_utils.make_dir(outfile)
     for of in outformat:
+        print 'Saving as {0}'.format(outfile+'.'+of)
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
 
 
-def myround(x, base=5, up=False, down=False):
-    if up == down and up is True: assert 0
-    if up: return int(base * np.round(float(x)/base-0.5))
-    elif down: return int(base * np.round(float(x)/base+0.5))
-    else: int(base * np.round(float(x)/base))
+def plot_sens_full(data, outfile, outformat, args):
+    print 'Making FULL sensitivity plot'
 
-
-def plot_BSM_angles_limit(dirname, outfile, outformat, args, bayesian):
-    """Make BSM angles vs scale limit plot."""
-    if not args.plot_angles_limit: return
-    print "Making BSM angles limit plot."""
-    fig_text = gen_figtext(args)
-    xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
-
-    raw = []
-    for root, dirs, filenames in os.walk(dirname):
-        for fn in filenames:
-            if fn[-4:] == '.npy':
-                raw.append(np.load(os.path.join(root, fn)))
-    raw = np.vstack(raw)
-    raw = raw[np.argsort(raw[:,0])]
-    print 'raw', raw
-    print 'raw.shape', raw.shape
-    sc_ranges = (
-        myround(np.min(raw[0][:,0]), up=True),
-        myround(np.max(raw[0][:,0]), down=True)
-    )
-
-    proc = []
-    if bayesian:
-        for idx, theta in enumerate(raw):
-            scale, evidences = theta.T
-            null = evidences[0]
-            reduced_ev = -(evidences - null)
-            al = scale[reduced_ev > np.log(10**(1/2.))]
-            proc.append((idx+1, al[0]))
-    else:
-        for idx, theta in enumerate(raw):
-            scale, llh = theta.T
-            delta_llh = -2 * (llh - np.max(llh))
-            # 90% CL for 1 dof
-            al = scale[delta_llh > 2.71]
-            proc.append((idx+1, al[0]))
-
-    limits = np.array(proc)
-    print 'limits',  limits
+    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
+    xticks = ['{0}'.format(x) for x in range(np.min(args.dimensions),
+                                             np.max(args.dimensions)+1)]
+    yranges = [np.inf, -np.inf]
+    legend_handles = []
 
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
 
-    ax.set_xticklabels(['']+xticks+[''])
-    ax.set_xlim(0, len(xticks)+1)
-    ax.set_ylim(sc_ranges[0], sc_ranges[-1])
-    ax.set_xlabel(r'BSM angle')
-    ylabel = r'${\rm log}_{10} \Lambda' + get_units(args.dimension) + r'$'
-    ax.set_ylabel(ylabel)
+    ax.set_xlim(args.dimensions[0]-1, args.dimensions[-1]+1)
+    ax.set_xticklabels([''] + xticks + [''])
+    ax.set_xlabel(r'BSM operator dimension ' + r'$d$')
+    ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$')
+
+    argsc = deepcopy(args)
+    for idim in xrange(len(data)):
+        dim = args.dimensions[idim]
+        print 'dim', dim
+        argsc.dimension = dim
+        for isrc in xrange(len(data[idim])):
+            src = args.source_ratios[isrc]
+            argsc.source_ratio = src
+
+            # fig_text = gen_figtext(argsc)
+            scales, statistic = data[idim][isrc].T
+            min_idx = np.argmin(scales)
+            null = statistic[min_idx]
+            if args.stat_method is StatCateg.BAYESIAN:
+                reduced_ev = -(statistic - null)
+                # al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief
+                al = scales[reduced_ev > 0.4] # Testing
+            elif args.stat_method is StatCateg.FREQUENTIST:
+                reduced_ev = -2*(statistic - null)
+                al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks
+            if len(al) == 0:
+                print 'No points for DIM {0} FRS {1} NULL {2}!'.format(
+                    dim, misc_utils.solve_ratio(src), null
+                )
+                print 'Reduced EV {0}'.format(reduced_ev)
+                continue
+            lim = al[0]
+            label = '[{0}, {1}, {2}]'.format(*misc_utils.solve_ratio(src))
+            if lim < yranges[0]: yranges[0] = lim
+            if lim > yranges[1]: yranges[1] = lim+4
+            line = plt.Line2D(
+                (dim-0.1, dim+0.1), (lim, lim), lw=3, color=colour[isrc], label=label
+            )
+            ax.add_line(line)
+            if idim == 0: legend_handles.append(line)
+            x_offset = isrc*0.05 - 0.05
+            ax.annotate(
+                s='', xy=(dim+x_offset, lim), xytext=(dim+x_offset, lim+3),
+                arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[isrc]}
+            )
 
-    for l in limits:
-        line = plt.Line2D(
-            (l[0]-0.1, l[0]+0.1), (l[1], l[1]), lw=3, color='r'
-        )
-        ax.add_line(line)
-        # ax.arrow(
-        #     l[0], l[1], 0, -1.5, head_width=0.05, head_length=0.2, fc='r',
-        #     ec='r', lw=2
-        # )
-        ax.annotate(
-            s='', xy=l, xytext=(l[0], l[1]+1.5),
-            arrowprops={'arrowstyle': '<-', 'lw': 1.5, 'color':'r'}
-        )
+    try:
+        yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
+        ax.set_ylim(yranges)
+    except: pass
 
+    ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right')
     for ymaj in ax.yaxis.get_majorticklocs():
-        ax.axhline(y=ymaj, ls='-', color='gray', alpha=0.4, linewidth=1)
+        ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1)
     for xmaj in ax.xaxis.get_majorticklocs():
-        ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1)
+        ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
 
+    misc_utils.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_sens_fixed_angle(data, outfile, outformat, args):
+    print 'Making FIXED_ANGLE sensitivity plot'
+
+    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
+    xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
+    argsc = deepcopy(args)
+    for idim in xrange(len(data)):
+        dim = args.dimensions[idim]
+        print '= dim', dim
+        argsc.dimension = dim
+
+        yranges = [np.inf, -np.inf]
+        legend_handles = []
+
+        fig = plt.figure(figsize=(7, 5))
+        ax = fig.add_subplot(111)
+        ax.set_xlim(0, len(xticks)+1)
+        ax.set_xticklabels([''] + xticks + [''])
+        ax.set_xlabel(r'BSM operator angle')
+        ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$')
+
+        for isrc in xrange(len(data[idim])):
+            src = args.source_ratios[isrc]
+            argsc.source_ratio = src
+            print '== src', src
+
+            for ian in xrange(len(data[idim][isrc])):
+                print '=== an', ian
+                scales, statistic = data[idim][isrc][ian].T
+                min_idx = np.argmin(scales)
+                null = statistic[min_idx]
+                if args.stat_method is StatCateg.BAYESIAN:
+                    reduced_ev = -(statistic - null)
+                    al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief
+                elif args.stat_method is StatCateg.FREQUENTIST:
+                    reduced_ev = -2*(statistic - null)
+                    al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks
+                print 'reduced_ev', reduced_ev
+                if len(al) == 0:
+                    print 'No points for DIM {0} FRS {1}\nreduced_ev {2}!'.format(
+                        dim, src, reduced_ev
+                    )
+                    continue
+                lim = al[0]
+                print 'limit = {0}'.format(lim)
+                label = '[{0}, {1}, {2}]'.format(*misc_utils.solve_ratio(src))
+                if lim < yranges[0]: yranges[0] = lim
+                if lim > yranges[1]: yranges[1] = lim+4
+                line = plt.Line2D(
+                    (ian+1-0.1, ian+1+0.1), (lim, lim), lw=3, color=colour[isrc], label=label
+                )
+                ax.add_line(line)
+                if len(legend_handles) < isrc+1:
+                    legend_handles.append(line)
+                x_offset = isrc*0.05 - 0.05
+                ax.annotate(
+                    s='', xy=(ian+1+x_offset, lim), xytext=(ian+1+x_offset, lim+3),
+                    arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[isrc]}
+                )
+
+        try:
+            yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
+            ax.set_ylim(yranges)
+        except: pass
+
+        ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right')
+        for ymaj in ax.yaxis.get_majorticklocs():
+            ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+        for xmaj in ax.xaxis.get_majorticklocs():
+            ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+
+        out = outfile + '_DIM{0}'.format(dim)
+        misc_utils.make_dir(out)
+        for of in outformat:
+            print 'Saving plot as {0}'.format(out+'.'+of)
+            fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150)
+
+
+        for ymaj in ax.yaxis.get_majorticklocs():
+            ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+        for xmaj in ax.xaxis.get_majorticklocs():
+            ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+
+
+def plot_sens_corr_angle(data, outfile, outformat, args):
+    print 'Making CORR_ANGLE sensitivity plot'
+
+    labels = [r'$sin^2(2\mathcal{O}_{12})$',
+              r'$sin^2(2\mathcal{O}_{13})$',
+              r'$sin^2(2\mathcal{O}_{23})$']
+
+    argsc = deepcopy(args)
+    for idim in xrange(len(data)):
+        dim = args.dimensions[idim]
+        print '= dim', dim
+        argsc.dimension = dim
+        for isrc in xrange(len(data[idim])):
+            src = args.source_ratios[isrc]
+            argsc.source_ratio = src
+            print '== src', src
+            for ian in xrange(len(data[idim][isrc])):
+                print '=== an', ian
+
+                d = data[idim][isrc][ian].reshape(
+                    len(data[idim][isrc][ian])**2, 3
+                )
+
+                fig = plt.figure(figsize=(7, 5))
+                ax = fig.add_subplot(111)
+                ax.set_ylim(0, 1)
+                ax.set_ylabel(labels[ian])
+                ax.set_xlabel(r'${\rm log}_{10} \Lambda^{-1}'+get_units(dim)+r'$')
+
+                xranges = [np.inf, -np.inf]
+                legend_handles = []
+
+                y = d[:,0]
+                x = d[:,1]
+                z = d[:,2]
+
+                print 'x', x
+                print 'y', y
+                print 'z', z
+                null_idx = np.argmin(x)
+                null = z[null_idx]
+                print 'null = {0}, for scale = {1}'.format(null, x[null_idx])
+
+                if args.stat_method is StatCateg.BAYESIAN:
+                    z_r = -(z - null)
+                elif args.stat_method is StatCateg.FREQUENTIST:
+                    z_r = -2*(z - null)
+                print 'scale', x
+                print 'reduced ev', z_r
+
+                pdat = np.array([x, y, z_r, np.ones(x.shape)]).T
+                print 'pdat', pdat
+                pdat_clean = []
+                for d in pdat:
+                    if not np.any(np.isnan(d)): pdat_clean.append(d)
+                pdat = np.vstack(pdat_clean)
+                sort_column = 3
+                pdat_sorted = pdat[pdat[:,sort_column].argsort()]
+                uni, c = np.unique(pdat[:,sort_column], return_counts=True)
+                print uni, c
+                print len(uni)
+                print np.unique(c)
+                
+                n = len(uni)
+                assert len(np.unique(c)) == 1
+                c = c[0]
+                col_array = []
+                for col in pdat_sorted.T:
+                    col_array.append(col.reshape(n, c))
+                col_array = np.stack(col_array)
+                sep_arrays = []
+                for x_i in xrange(n):
+                    sep_arrays.append(col_array[:,x_i])
+                
+                print len(sep_arrays)
+                sep_arrays = sep_arrays[0][:3]
+                print sep_arrays
+
+                if args.stat_method is StatCateg.BAYESIAN:
+                    reduced_pdat_mask = (sep_arrays[2] > log(10**(3/2.))) # Strong degree of belief
+                elif args.stat_method is StatCateg.FREQUENTIST:
+                    reduced_pdat_mask = (sep_arrays[2] > 4.61) # 90% CL for 2 DOFS via Wilks
+                reduced_pdat = sep_arrays.T[reduced_pdat_mask].T
+                print 'reduced_pdat', reduced_pdat
+
+                ax.tick_params(axis='x', labelsize=11)
+                ax.tick_params(axis='y', labelsize=11)
+
+                mini, maxi = np.min(x), np.max(x)
+                ax.set_xlim((mini, maxi))
+                ax.set_ylim(0, 1)
+                ax.grid(b=False)
+                
+                x_v = reduced_pdat[0].round(decimals=4)
+                y_v = reduced_pdat[1].round(decimals=4)
+                uniques = np.unique(x_v)
+                print 'uniques', uniques
+                if len(uniques) == 1: continue
+                bw = np.min(np.diff(uniques))
+                print 'bw', bw
+                print np.diff(uniques)
+                uni_x_split = np.split(uniques, np.where(np.diff(uniques) > bw*(1e20))[0] + 1)
+                print 'len(uni_x_split)', len(uni_x_split)
+                for uni_x in uni_x_split:
+                    p_x_l, p_y_l = [], []
+                    p_x_u, p_y_u = [], []
+                    for uni in uni_x:
+                        idxes = np.where(x_v == uni)[0]
+                        ymin, ymax = 1, 0
+                        for idx in idxes:
+                            if y_v[idx] < ymin: ymin = y_v[idx]
+                            if y_v[idx] > ymax: ymax = y_v[idx]
+                        p_x_l.append(uni)
+                        p_y_l.append(ymin)
+                        p_x_u.append(uni)
+                        p_y_u.append(ymax)
+                    p_x_l, p_y_l = np.array(p_x_l, dtype=np.float64), np.array(p_y_l, dtype=np.float64)
+                    p_x_u, p_y_u = np.array(list(reversed(p_x_u)), dtype=np.float64), np.array(list(reversed(p_y_u)), dtype=np.float64)
+                    p_x = np.hstack([p_x_l, p_x_u])
+                    p_y = np.hstack([p_y_l, p_y_u])
+                    p_x = np.r_[p_x, p_x[0]]
+                    p_y = np.r_[p_y, p_y[0]]
+                    print 'p_x', p_x
+                    print 'p_y', p_y
+                    try:
+                        tck, u = interpolate.splprep([p_x, p_y], s=1e-5, per=True)
+                        xi, yi = interpolate.splev(np.linspace(0, 1, 1000), tck)
+                    except:
+                        xi, yi = p_x, p_y
+                    ax.fill(xi, yi, 'r', edgecolor='r', linewidth=1)
+
+                out = outfile + '_DIM{0}_SRC{1}_AN{2}'.format(dim, isrc, ian)
+                misc_utils.make_dir(out)
+                for of in outformat:
+                    print 'Saving plot as {0}'.format(out+'.'+of)
+                    fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150)