Wed May 23 16:23:12 CDT 2018

1 files changed, 0 insertions, 414 deletions

diff --git a/plot_bf.py b/plot_bf.py
deleted file mode 100755
index 13664b9..0000000
--- a/plot_bf.py
+++ /dev/null
@@ -1,414 +0,0 @@
-#! /usr/bin/env python
-# author : S. Mandalia
-#          s.p.mandalia@qmul.ac.uk
-#
-# date   : April 14, 2018
-
-"""
-HESE BSM flavour ratio sensivity plotting script
-"""
-
-from __future__ import absolute_import, division
-
-import os
-
-import numpy as np
-import numpy.ma as ma
-
-import matplotlib as mpl
-mpl.use('Agg')
-from matplotlib import rc
-from matplotlib import pyplot as plt
-from matplotlib.offsetbox import AnchoredText
-
-import fr
-from utils import misc as misc_utils
-from utils.fr import normalise_fr
-from utils.plot import bayes_factor_plot, myround, get_units
-
-
-rc('text', usetex=False)
-rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
-
-fix_sfr_mfr = [
-    (1, 1, 1, 1, 2, 0),
-    # (1, 1, 1, 1, 0, 0),
-    (1, 1, 1, 0, 1, 0),
-]
-
-# FR
-dimension         = [3, 6]
-# dimension         = [3, 4, 5, 6, 7, 8]
-sigma_ratio       = ['0.01']
-energy_dependance = 'spectral'
-spectral_index    = -2
-binning           = [1e4, 1e7, 5]
-fix_mixing        = 'False'
-fix_mixing_almost = 'False'
-scale_region      = "1E10"
-
-# Likelihood
-likelihood = 'golemfit'
-
-# Nuisance
-convNorm        = 1.
-promptNorm      = 0.
-muonNorm        = 1.
-astroNorm       = 6.9
-astroDeltaGamma = 2.5
-
-# GolemFit
-ast  = 'p2_0'
-data = 'real'
-
-# Bayes Factor
-bayes_bins        = 100
-bayes_live_points = 1000
-bayes_tolerance   = 0.01
-bayes_eval_bin    = True # set to 'all' to run normally
-
-# Plot
-plot_bayes        = False
-plot_angles_limit = True
-plot_angles_corr  = False
-outformat         = ['png']
-# significance      = np.log(10**(1/2.))
-significance      = np.log(10**(3/2.))
-
-
-bayes_array = ma.masked_equal(np.zeros((len(dimension), len(fix_sfr_mfr), bayes_bins, 2)), 0)
-angles_lim_array = np.zeros((len(dimension), len(fix_sfr_mfr), 3, bayes_bins, 2))
-angles_corr_array = np.zeros((len(dimension), len(fix_sfr_mfr), 3, bayes_bins, bayes_bins, 3))
-for i_dim, dim in enumerate(dimension):
-    if energy_dependance == 'mono':
-        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/{2:.0E}'.format(likelihood, dim, en)
-    elif energy_dependance == 'spectral':
-        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(likelihood, dim, spectral_index)
-
-    bayes_output = 'None'
-    angles_lim_output = 'None'
-    angles_corr_output = 'None'
-    for sig in sigma_ratio:
-        for i_frs, frs in enumerate(fix_sfr_mfr):
-            outchains = outchain_head + '/fix_ifr/{0}/'.format(str(sig).replace('.', '_'))
-            if plot_bayes:
-                bayes_output = outchains + '/bayes_factor/'
-            if plot_angles_limit:
-                angles_lim_output = outchains + '/angles_limit/'
-            if plot_angles_corr:
-                angles_corr_output = outchains + '/angles_corr/'
-
-            argstring = '--measured-ratio {0} {1} {2} --fix-source-ratio True --source-ratio {3} {4} {5} --dimension {6} --seed 24 --outfile {7} --run-mcmc False --likelihood {8} --plot-angles False --bayes-output {9} --angles-lim-output {10} --bayes-bins {11} --angles-corr-output {12}'.format(frs[0], frs[1], frs[2], frs[3], frs[4], frs[5], dim, outchains, likelihood, bayes_output, angles_lim_output, bayes_bins, angles_corr_output)
-            args = fr.parse_args(argstring)
-            fr.process_args(args)
-            # misc_utils.print_args(args)
-
-            if plot_bayes:
-                infile = args.bayes_output+'/fr_evidence'+misc_utils.gen_identifier(args)
-            if plot_angles_limit:
-                infile = args.angles_lim_output+'/fr_an_evidence'+misc_utils.gen_identifier(args)
-            if plot_angles_corr:
-                infile = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args)
-            scan_scales = np.linspace(
-                np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.bayes_bins
-            )
-            print 'scan_scales', scan_scales
-            raw = []
-            fail = 0
-            if plot_angles_corr:
-                scan_angles = np.linspace(0, 1, args.bayes_bins)
-            for i_sc, sc in enumerate(scan_scales):
-                if plot_angles_corr:
-                    for i_an, an in enumerate(scan_angles):
-                        idx = i_sc*args.bayes_bins + i_an
-                        infile_s = infile + '_idx_{0}'.format(idx)
-                        try:
-                            lf = np.load(infile_s+'.npy')
-                            print 'lf.shape', lf.shape
-                        except IOError:
-                            fail += 1
-                            print 'failed to open {0}'.format(infile_s)
-                            lf = np.full((3, 1, 1, 3), np.nan)
-                            pass
-                        for x in xrange(len(lf)):
-                            angles_corr_array[i_dim][i_frs][x][i_sc][i_an] = np.array(lf[x])
-                    continue
-                try:
-                    infile_s = infile + '_scale_{0:.0E}'.format(np.power(10, sc))
-                    lf = np.load(infile_s+'.npy')
-                    print lf.shape
-                    if plot_angles_limit:
-                        if len(lf.shape) == 3: lf = lf[:,0,:]
-                    raw.append(lf)
-                except IOError:
-                    fail += 1
-                    print 'failed to open {0}'.format(infile_s)
-                    if plot_bayes:
-                        raw.append([0, 0])
-                    if plot_angles_limit:
-                        raw.append(np.zeros((3, 2)))
-                    pass
-            print 'failed to open {0} files'.format(fail)
-
-            if plot_bayes:
-                raw = np.vstack(raw)
-            if plot_angles_limit:
-                raw = np.vstack(raw).reshape(args.bayes_bins, 3, 2)
-                a = ma.masked_equal(np.zeros((3, args.bayes_bins, 2)), 0)
-                for i_x, x in enumerate(raw):
-                    for i_y, y in enumerate(x):
-                        a[i_y][i_x] = ma.masked_equal(y, 0)
-            if plot_angles_corr:
-                a = angles_corr_array[i_dim][i_frs]
-                a = ma.masked_invalid(a, 0)
-                # for i_sc in xrange(len(scan_scales)):
-                #     for i_a in xrange(len(scan_angles)):
-                #         try:
-                #             bayes_factor_plot(
-                #                 a[i_sc,:,i_a,:][:,(0,2)], './mnrun/corr/test_corr_DIM{0}_FR{1}_AN{2}_SC{3}'.format(dim, i_frs, i_a, i_sc), ['png'], args
-                #             )
-                #         except: pass
-
-            if plot_bayes:
-                bayes_array[i_dim][i_frs] = ma.masked_equal(raw, 0)
-                bayes_factor_plot( bayes_array[i_dim][i_frs], './mnrun/test_full_DIM{0}_FR{1}'.format(dim, i_frs), ['png'], args
-                )
-
-            if plot_angles_limit:
-                angles_lim_array[i_dim][i_frs] = ma.masked_equal(a, 0)
-                for i_a, angle in enumerate(a):
-                    bayes_factor_plot(
-                        angle, './mnrun/test_angles_DIM{0}_FR{1}_AN{2}'.format(i_dim, i_frs, i_a), ['png'], args
-                    )
-
-if plot_bayes:
-    fig = plt.figure(figsize=(7, 5))
-    ax = fig.add_subplot(111)
-
-    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-    yranges = [np.inf, -np.inf]
-    legend_handles = []
-    ax.set_xlim(dimension[0]-1, dimension[-1]+1)
-    xticks = [''] + range(dimension[0], dimension[-1]+1) + ['']
-    ax.set_xticklabels(xticks)
-    ax.set_xlabel(r'BSM operator dimension ' + r'$d$')
-    ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$')
-    for i_dim, dim in enumerate(dimension):
-        for i_frs, frs in enumerate(fix_sfr_mfr):
-            scale, evidences = bayes_array[i_dim][i_frs].T
-            null = evidences[np.argmin(scale)]
-            # TODO(shivesh): negative or not?
-            reduced_ev = -(evidences - null)
-            al = scale[reduced_ev > significance]
-            if len(al) > 0:
-                label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5])
-                lim = al[0]
-                print 'frs, dim, lim = ', frs, dim, lim
-                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[i_frs], label=label
-                )
-                ax.add_line(line)
-                if i_dim == 0: legend_handles.append(line)
-                x_offset = i_frs*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[i_frs]}
-                )
-
-            else:
-                print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, null)
-                # print 'scales, reduced_ev', np.dstack([scale.data, reduced_ev.data])
-    try:
-        yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
-        ax.set_ylim(yranges)
-    except: pass
-
-    ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right')
-    for ymaj in ax.yaxis.get_majorticklocs():
-        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)
-
-    for of in outformat:
-        fig.savefig('./images/bayes/bayes_factor.'+of, bbox_inches='tight', dpi=150)
-
-if plot_angles_limit:
-    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-    for i_dim, dim in enumerate(dimension):
-        fig = plt.figure(figsize=(7, 5))
-        ax = fig.add_subplot(111)
-        yranges = [np.inf, -np.inf]
-        legend_handles = []
-        xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
-        ax.set_xlim(0, len(xticks)+1)
-        ax.set_xticklabels([''] + xticks + [''])
-        ax.set_xlabel(r'BSM operator angle')
-        ylabel = r'${\rm log}_{10} \Lambda^{-1}' + get_units(dim) + r'$'
-        ax.set_ylabel(ylabel)
-        for i_th in xrange(len(xticks)):
-            for i_frs, frs in enumerate(fix_sfr_mfr):
-                scale, evidences = angles_lim_array[i_dim][i_frs][i_th].T
-                null = evidences[np.argmin(scale)]
-                # TODO(shivesh): negative or not?
-                reduced_ev = -(evidences - null)
-                al = scale[reduced_ev > significance]
-                # print 'scales, reduced_ev', np.dstack([scale, reduced_ev])
-                if len(al) > 0:
-                    label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5])
-                    lim = al[0]
-                    print 'frs, dim, lim = ', frs, dim, lim
-                    if lim < yranges[0]: yranges[0] = lim
-                    if lim > yranges[1]: yranges[1] = lim+4
-                    line = plt.Line2D(
-                        (i_th+1-0.1, i_th+1+0.1), (lim, lim), lw=3, color=colour[i_frs], label=label
-                    )
-                    ax.add_line(line)
-                    if i_th == 0: legend_handles.append(line)
-                    x_offset = i_frs*0.05 - 0.05
-                    ax.annotate(
-                        s='', xy=(i_th+1+x_offset, lim), xytext=(i_th+1+x_offset, lim+3),
-                        arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[i_frs]}
-                    )
-        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',
-                  title='dimension {0}'.format(dim))
-
-        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)
-
-        for of in outformat:
-            fig.savefig('./images/bayes/angles_limit_DIM{0}'.format(dim)+'.'+of, bbox_inches='tight', dpi=150)
-
-if plot_angles_corr:
-    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-    labels = [r'$sin^2(2\mathcal{O}_{12})$', r'$sin^2(2\mathcal{O}_{13})$', r'$sin^2(2\mathcal{O}_{23})$']
-    for i_dim, dim in enumerate(dimension):
-        for i_frs, frs in enumerate(fix_sfr_mfr):
-            print '== DIM{0}'.format(dim)
-            print '== FRS = {0}'.format(frs)
-            array = angles_corr_array[i_dim][i_frs]
-            print 'array', array
-            print 'array.shape', array.shape
-            for i_scen in xrange(len(labels)):
-                d = array[i_scen].reshape(len(array[i_scen])**2, 3)
-                print 'd.shape', d.shape
-                fig = plt.figure(figsize=(7, 5))
-                ax = fig.add_subplot(111)
-                xranges = [np.inf, -np.inf]
-                legend_handles = []
-                ax.set_ylim(0, 1)
-                ax.set_ylabel(labels[i_scen])
-                xlabel = r'${\rm log}_{10} \Lambda^{-1}' + get_units(dim) + r'$'
-                ax.set_xlabel(xlabel)
-
-                x = d[:,0]
-                y = d[:,1]
-                z = d[:,2]
-
-                data_clean = []
-                for id in d:
-                    if not np.any(np.isnan(id)): data_clean.append(id)
-                d_c = np.vstack(data_clean)
-
-                x = d_c[:,0]
-                y = d_c[:,1]
-                z = d_c[:,2]
-
-                # print 'x', x
-                # print 'y', y
-                # print 'z', z
-                null_idx = np.argmin(x)
-                null = z[null_idx]
-                print 'null = {0}, for scale = {1}'.format(null, x[null_idx])
-                z = -(z - null)
-                print 'scale', x
-                print 'bayes_factor', z
-
-                # x_ = np.linspace(np.min(x), np.max(x), 30)
-                # y_ = np.linspace(np.min(y), np.max(y), 30)
-                # z_ = interpolate.gridddata((x, y), z, (x_[None,:], y_[:,None]), method='nearest')
-
-                data = np.array([x, y, z, np.ones(x.shape)]).T
-                sort_column = 3
-                data_sorted = data[data[:,sort_column].argsort()]
-                uni, c = np.unique(data[:,sort_column], return_counts=True)
-                print uni, c
-                print len(uni)
-                print np.unique(c)
-                
-                n = len(uni)
-                assert len(np.unique(c)) == 1
-                c = c[0]
-                col_array = []
-                for col in data_sorted.T:
-                    col_array.append(col.reshape(n, c))
-                col_array = np.stack(col_array)
-                sep_arrays = []
-                for x_i in xrange(n):
-                    sep_arrays.append(col_array[:,x_i])
-                
-                print len(sep_arrays)
-                sep_arrays = sep_arrays[0][:3]
-                print sep_arrays
-
-                allowed_bf = (sep_arrays[2] < significance) # Shade the excluded region
-                data_allowed_bf = sep_arrays.T[allowed_bf].T
-                print 'data_allowed_bf', data_allowed_bf
-
-                ax.tick_params(axis='x', labelsize=11)
-                ax.tick_params(axis='y', labelsize=11)
-
-                mini, maxi = np.min(scan_scales), np.max(scan_scales)
-                ax.set_xlim((mini, maxi))
-                ax.set_ylim(0, 1)
-                ax.grid(b=False)
-                
-                x_v = data_allowed_bf[0].round(decimals=4)
-                y_v = data_allowed_bf[1].round(decimals=4)
-                uniques = np.unique(x_v)
-                # print 'uniques', uniques
-                if len(uniques) == 1: continue
-                bw = np.min(np.diff(uniques))
-                # print 'bw', bw
-                print np.diff(uniques)
-                uni_x_split = np.split(uniques, np.where(np.diff(uniques) > bw*(1e20))[0] + 1)
-                # print 'len(uni_x_split)', len(uni_x_split)
-                for uni_x in uni_x_split:
-                    p_x_l, p_y_l = [], []
-                    p_x_u, p_y_u = [], []
-                    for uni in uni_x:
-                        idxes = np.where(x_v == uni)[0]
-                        ymin, ymax = 1, 0
-                        for idx in idxes:
-                            if y_v[idx] < ymin: ymin = y_v[idx]
-                            if y_v[idx] > ymax: ymax = y_v[idx]
-                        p_x_l.append(uni)
-                        p_y_l.append(ymin)
-                        p_x_u.append(uni)
-                        p_y_u.append(ymax)
-                    p_x_l, p_y_l = np.array(p_x_l, dtype=np.float64), np.array(p_y_l, dtype=np.float64)
-                    p_x_u, p_y_u = np.array(list(reversed(p_x_u)), dtype=np.float64), np.array(list(reversed(p_y_u)), dtype=np.float64)
-                    p_x = np.hstack([p_x_l, p_x_u])
-                    p_y = np.hstack([p_y_l, p_y_u])
-                    p_x = np.r_[p_x, p_x[0]]
-                    p_y = np.r_[p_y, p_y[0]]
-                    print 'p_x', p_x
-                    print 'p_y', p_y
-                    try:
-                        tck, u = interpolate.splprep([p_x, p_y], s=1e-5, per=True)
-                        xi, yi = interpolate.splev(np.linspace(0, 1, 1000), tck)
-                    except:
-                        xi, yi = p_x, p_y
-                    ax.fill(xi, yi, 'r', edgecolor='r', linewidth=1)
-                
-                for of in outformat:
-                    plt.savefig('./images/bayes/lim_corr_DIM{0}_AN{1}_FRS{2}'.format(dim, i_scen, i_frs)+of, bbox_inches='tight', dpi=150)
-