diff options
| -rw-r--r-- | bout/plot.py | 129 | ||||
| -rw-r--r-- | bout/plot_corr.py | 193 | ||||
| -rw-r--r-- | bout/plot_full.py | 127 | ||||
| -rwxr-xr-x | fr.py | 481 | ||||
| -rwxr-xr-x | fr_edit.py | 581 | ||||
| -rwxr-xr-x[-rw-r--r--] | plot_bf.py | 204 | ||||
| -rwxr-xr-x | sens.py | 480 | ||||
| -rw-r--r-- | submitter/make_dag.py | 21 | ||||
| -rw-r--r-- | test/LV.out | 74 | ||||
| -rw-r--r-- | test/NSI.out | 75 | ||||
| -rw-r--r-- | test/test.png | bin | 104565 -> 0 bytes | |||
| -rw-r--r-- | test/test_LV.py | 95 | ||||
| -rw-r--r-- | test/test_NSI.png | bin | 117422 -> 0 bytes | |||
| -rw-r--r-- | test/test_NSI.py | 103 | ||||
| -rw-r--r-- | utils/enums.py | 11 | ||||
| -rw-r--r-- | utils/fr.py | 90 | ||||
| -rw-r--r-- | utils/gf.py | 1 | ||||
| -rw-r--r-- | utils/likelihood.py | 33 | ||||
| -rw-r--r-- | utils/mcmc.py | 9 | ||||
| -rw-r--r-- | utils/misc.py | 242 | ||||
| -rw-r--r-- | utils/multinest.py | 117 | ||||
| -rw-r--r-- | utils/param.py | 235 | ||||
| -rw-r--r-- | utils/plot.py | 112 |
23 files changed, 2267 insertions, 1146 deletions
diff --git a/bout/plot.py b/bout/plot.py new file mode 100644 index 0000000..717cf81 --- /dev/null +++ b/bout/plot.py @@ -0,0 +1,129 @@ +import os + +import numpy as np +import numpy.ma as ma + +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import pyplot as plt +from matplotlib.offsetbox import AnchoredText +from matplotlib import rc + +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, 6] +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' +confidence = 2.71 # 90% for 1DOF +outformat = ['png'] + + +def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01): + mr = np.array(measured_ratio) / float(np.sum(measured_ratio)) + sr = np.array(source_ratio) / float(np.sum(source_ratio)) + si = sigma_ratio + out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format( + int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), + int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension + ) + return out + + +def get_units(dimension): + if dimension == 3: return r' / GeV' + if dimension == 4: return r'' + if dimension == 5: return r' / GeV^{-1}' + if dimension == 6: return r' / GeV^{-2}' + if dimension == 7: return r' / GeV^{-3}' + if dimension == 8: return r' / GeV^{-4}' + + +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)) + + +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' + get_units(dim) + r'$' + ax.set_ylabel(ylabel) + for i_frs, frs in enumerate(fix_sfr_mfr): + print '== DIM{0}'.format(dim) + print '== FRS = {0}'.format(frs) + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index) + infile = outchain_head + '/angles_limit/fr_anfr_evidence'+ gen_identifier(frs[:3], frs[-3:], dim) + '.npy' + try: + array = np.load(infile) + except IOError: + print 'failed to open {0}'.format(infile) + continue + print 'array', array + print 'array', array.shape + for i_th in xrange(len(xticks)): + scale, llhs = array[i_th].T + min_llh = np.min(llhs) + delta_llh = 2*(llhs - min_llh) + print 'scale', scale + print 'delta_llh', delta_llh + al = scale[delta_llh < confidence] + 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]} + ) + else: + print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, min_llh) + try: + yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) + # ax.set_ylim(yranges) + ax.set_ylim([-30, -20]) + 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/freq/lim_DIM{0}.'.format(dim)+of, bbox_inches='tight', dpi=150) diff --git a/bout/plot_corr.py b/bout/plot_corr.py new file mode 100644 index 0000000..a75fe8a --- /dev/null +++ b/bout/plot_corr.py @@ -0,0 +1,193 @@ +import os + +import numpy as np +import numpy.ma as ma + +import scipy.interpolate as interpolate + +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import pyplot as plt +from matplotlib.offsetbox import AnchoredText +from matplotlib import rc + +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, 6] +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' +confidence = 4.61 # 90% for 2DOF +outformat = ['png'] + + +def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01): + mr = np.array(measured_ratio) / float(np.sum(measured_ratio)) + sr = np.array(source_ratio) / float(np.sum(source_ratio)) + si = sigma_ratio + out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format( + int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), + int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension + ) + return out + + +def get_units(dimension): + if dimension == 3: return r' / GeV' + if dimension == 4: return r'' + if dimension == 5: return r' / GeV^{-1}' + if dimension == 6: return r' / GeV^{-2}' + if dimension == 7: return r' / GeV^{-3}' + if dimension == 8: return r' / GeV^{-4}' + + +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)) + + +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) + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index) + infile = outchain_head + '/angles_corr/fr_co_evidence'+ gen_identifier(frs[:3], frs[-3:], dim) + '.npy' + # infile = '../mnrun/fr_co_evidence_033_033_033_0010_sfr_033_066_000_DIM6_single_scale.npy' + try: + array = ma.masked_invalid(np.load(infile)) + except IOError: + print 'failed to open {0}'.format(infile) + continue + print 'array', array + print 'array', array.shape + for i_scen in xrange(len(labels)): + d = array[i_scen].reshape(len(array[i_scen])**2, 3) + 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' + get_units(dim) + r'$' + ax.set_xlabel(xlabel) + + x = d[:,0] + y = d[:,1] + z = d[:,2] + + print 'x', x + print 'y', y + print 'z', z + null_idx = np.argmin(z) + null = z[null_idx] + print 'null = {0}, for scale = {1}'.format(null, x[null_idx]) + z = 2*(z - null) + print 'scale', x + print 'delta_llh', 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 + print 'data', data + data_clean = [] + for d in data: + if not np.any(np.isnan(d)): data_clean.append(d) + data = np.vstack(data_clean) + 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 + + llh_90_percent = (sep_arrays[2] < confidence) + data_90_percent = sep_arrays.T[llh_90_percent].T + print 'data_90_percent', data_90_percent + + 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 = data_90_percent[0].round(decimals=4) + y_v = data_90_percent[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/freq/lim_corr_DIM{0}_AN{1}_FRS{2}'.format(dim, i_scen, i_frs)+of, bbox_inches='tight', dpi=150) + diff --git a/bout/plot_full.py b/bout/plot_full.py new file mode 100644 index 0000000..f2e1919 --- /dev/null +++ b/bout/plot_full.py @@ -0,0 +1,127 @@ +import os + +import numpy as np +import numpy.ma as ma + +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import pyplot as plt +from matplotlib.offsetbox import AnchoredText +from matplotlib import rc + +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' +confidence = 2.71 # 90% for 1DOF +outformat = ['png'] + + +def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01): + mr = np.array(measured_ratio) / float(np.sum(measured_ratio)) + sr = np.array(source_ratio) / float(np.sum(source_ratio)) + si = sigma_ratio + out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format( + int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), + int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension + ) + return out + + +def get_units(dimension): + if dimension == 3: return r' / GeV' + if dimension == 4: return r'' + if dimension == 5: return r' / GeV^{-1}' + if dimension == 6: return r' / GeV^{-2}' + if dimension == 7: return r' / GeV^{-3}' + if dimension == 8: return r' / GeV^{-4}' + + +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)) + + +colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + +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 / GeV^{-d+4}$') +for i_dim, dim in enumerate(dimension): + for i_frs, frs in enumerate(fix_sfr_mfr): + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index) + infile = outchain_head + '/bayes_factor/fr_fr_evidence' + gen_identifier(frs[:3], frs[-3:], dim) + '.npy' + try: + array = np.load(infile) + except IOError: + print 'failed to open {0}'.format(infile) + continue + print 'array', array + print 'array', array.shape + scale, llhs = array.T + print 'scale min', scale[np.argmin(llhs)] + null = llhs[np.argmin(llhs)] + # null = llhs[0] + # TODO(shivesh): negative or not? + reduced_ev = 2*(llhs - null) + print 'reduced_ev', reduced_ev + al = scale[reduced_ev < confidence] + 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]) +yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) +ax.set_ylim(yranges) + +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/freq/full_corr.'+of, bbox_inches='tight', dpi=150) @@ -5,27 +5,25 @@ # date : March 17, 2018 """ -HESE BSM flavour ratio analysis script +HESE BSM flavour ratio MCMC analysis script """ from __future__ import absolute_import, division -import os import argparse from functools import partial import numpy as np -import GolemFitPy as gf - +from utils import fr as fr_utils from utils import gf as gf_utils from utils import likelihood as llh_utils from utils import mcmc as mcmc_utils from utils import misc as misc_utils from utils import plot as plot_utils from utils.enums import EnergyDependance, Likelihood, MCMCSeedType, ParamTag -from utils.fr import MASS_EIGENVALUES, normalise_fr, fr_to_angles -from utils.misc import enum_parse, Param, ParamSet +from utils.fr import estimate_scale, normalise_fr +from utils.param import Param, ParamSet, get_paramsets def define_nuisance(): @@ -52,56 +50,6 @@ def nuisance_argparse(parser): ) -def get_paramsets(args): - """Make the paramsets for generating the Asmimov MC sample and also running - the MCMC. - """ - asimov_paramset = [] - mcmc_paramset = [] - if args.likelihood == Likelihood.GOLEMFIT: - nuisance_paramset = define_nuisance() - asimov_paramset.extend(nuisance_paramset.params) - mcmc_paramset.extend(nuisance_paramset.params) - for parm in nuisance_paramset: - parm.value = args.__getattribute__(parm.name) - tag = ParamTag.BESTFIT - flavour_angles = fr_to_angles(args.measured_ratio) - asimov_paramset.extend([ - Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag), - Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag), - ]) - asimov_paramset = ParamSet(asimov_paramset) - - if not args.fix_mixing and not args.fix_mixing_almost: - tag = ParamTag.MMANGLES - e = 1e-10 - mcmc_paramset.extend([ - Param(name='s_12^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag), - Param(name='c_13^4', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag), - Param(name='s_23^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag), - Param(name='dcp', value=np.pi, ranges=[0.+e, 2*np.pi-e], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag) - ]) - if args.fix_mixing_almost: - tag = ParamTag.MMANGLES - mcmc_paramset.extend([ - Param(name='s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag) - ]) - if not args.fix_scale: - tag = ParamTag.SCALE - mcmc_paramset.append( - Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\Lambda'+plot_utils.get_units(args.dimension), tag=tag) - ) - if not args.fix_source_ratio: - tag = ParamTag.SRCANGLES - mcmc_paramset.extend([ - Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag), - Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag) - ]) - mcmc_paramset = ParamSet(mcmc_paramset) - return asimov_paramset, mcmc_paramset - - def process_args(args): """Process the input args.""" if args.fix_mixing and args.fix_scale: @@ -110,10 +58,6 @@ def process_args(args): raise NotImplementedError( '--fix-mixing and --fix-mixing-almost cannot be used together' ) - if args.run_bayes_factor and args.fix_scale: - raise NotImplementedError( - '--run-bayes-factor and --fix-scale cannot be used together' - ) args.measured_ratio = normalise_fr(args.measured_ratio) if args.fix_source_ratio: @@ -125,26 +69,9 @@ def process_args(args): ) if not args.fix_scale: - if args.energy_dependance is EnergyDependance.MONO: - args.scale = np.power( - 10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \ - np.log10(args.energy**(args.dimension-3)) + 6 - ) - elif args.energy_dependance is EnergyDependance.SPECTRAL: - args.scale = np.power( - 10, np.round( - np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \ - - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3)) - ) + 6 - ) - """Estimate the scale of when to expect the BSM term to contribute.""" + args.scale = estimate_scale(args) args.scale_region = (args.scale/args.scale_region, args.scale*args.scale_region) - if args.bayes_eval_bin.lower() == 'all': - args.bayes_eval_bin = None - else: - args.bayes_eval_bin = int(args.bayes_eval_bin) - def parse_args(args=None): """Parse command line arguments""" @@ -153,108 +80,7 @@ def parse_args(args=None): formatter_class=misc_utils.SortingHelpFormatter, ) parser.add_argument( - '--measured-ratio', type=int, nargs=3, default=[1, 1, 1], - help='Set the central value for the measured flavour ratio at IceCube' - ) - parser.add_argument( - '--sigma-ratio', type=float, default=0.01, - help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH' - ) - parser.add_argument( - '--fix-source-ratio', type=misc_utils.parse_bool, default='False', - help='Fix the source flavour ratio' - ) - parser.add_argument( - '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance), - choices=EnergyDependance, - help='Type of energy dependance to use in the BSM fit' - ) - parser.add_argument( - '--spectral-index', default=-2, type=int, - help='Spectral index for spectral energy dependance' - ) - parser.add_argument( - '--binning', default=[1e4, 1e7, 5], type=float, nargs=3, - help='Binning for spectral energy dependance' - ) - parser.add_argument( - '--run-bayes-factor', type=misc_utils.parse_bool, default='False', - help='Make the bayes factor plot for the scale' - ) - parser.add_argument( - '--run-angles-limit', type=misc_utils.parse_bool, default='False', - help='Make the limit vs BSM angles plot' - ) - parser.add_argument( - '--run-angles-correlation', type=misc_utils.parse_bool, default='False', - help='Make the limit vs BSM angles plot' - ) - parser.add_argument( - '--bayes-bins', type=int, default=10, - help='Number of bins for the Bayes factor plot' - ) - parser.add_argument( - '--bayes-eval-bin', type=str, default='all', - help='Which bin to evalaute for Bayes factor plot' - ) - parser.add_argument( - '--bayes-live-points', type=int, default=400, - help='Number of live points for MultiNest evaluations' - ) - parser.add_argument( - '--bayes-tolerance', type=float, default=0.5, - help='Tolerance for MultiNest' - ) - parser.add_argument( - '--bayes-output', type=str, default='./mnrun/', - help='Folder to store MultiNest evaluations' - ) - parser.add_argument( - '--angles-lim-output', type=str, default='./mnrun/', - help='Folder to store MultiNest evaluations' - ) - parser.add_argument( - '--angles-corr-output', type=str, default='./mnrun/', - help='Folder to store MultiNest evaluations' - ) - parser.add_argument( - '--source-ratio', type=int, nargs=3, default=[2, 1, 0], - help='Set the source flavour ratio for the case when you want to fix it' - ) - parser.add_argument( - '--no-bsm', type=misc_utils.parse_bool, default='False', - help='Turn off BSM terms' - ) - parser.add_argument( - '--fix-mixing', type=misc_utils.parse_bool, default='False', - help='Fix all mixing parameters to bi-maximal mixing' - ) - parser.add_argument( - '--fix-mixing-almost', type=misc_utils.parse_bool, default='False', - help='Fix all mixing parameters except s23' - ) - parser.add_argument( - '--fix-scale', type=misc_utils.parse_bool, default='False', - help='Fix the new physics scale' - ) - parser.add_argument( - '--scale', type=float, default=1e-30, - help='Set the new physics scale' - ) - parser.add_argument( - '--scale-region', type=float, default=1e7, - help='Set the size of the box to scan for the scale' - ) - parser.add_argument( - '--dimension', type=int, default=3, - help='Set the new physics dimension to consider' - ) - parser.add_argument( - '--energy', type=float, default=1000, - help='Set the energy scale' - ) - parser.add_argument( - '--seed', type=int, default=99, + '--seed', type=misc_utils.seed_parse, default='25', help='Set the random seed value' ) parser.add_argument( @@ -263,12 +89,12 @@ def parse_args(args=None): ) parser.add_argument( '--outfile', type=str, default='./untitled', - help='Path to output chains' + help='Path to output results' ) + fr_utils.fr_argparse(parser) + gf_utils.gf_argparse(parser) llh_utils.likelihood_argparse(parser) mcmc_utils.mcmc_argparse(parser) - gf_utils.gf_argparse(parser) - plot_utils.plot_argparse(parser) nuisance_argparse(parser) if args is None: return parser.parse_args() else: return parser.parse_args(args.split()) @@ -279,31 +105,31 @@ def main(): process_args(args) misc_utils.print_args(args) - np.random.seed(args.seed) + if args.seed is not None: + np.random.seed(args.seed) - asimov_paramset, mcmc_paramset = get_paramsets(args) + asimov_paramset, llh_paramset = get_paramsets(args) outfile = misc_utils.gen_outfile_name(args) print '== {0:<25} = {1}'.format('outfile', outfile) if args.run_mcmc: if args.likelihood is Likelihood.GOLEMFIT: fitter = gf_utils.setup_fitter(args, asimov_paramset) - else: - fitter = None + else: fitter = None ln_prob = partial( llh_utils.ln_prob, args=args, fitter=fitter, - asimov_paramset=asimov_paramset, mcmc_paramset=mcmc_paramset + asimov_paramset=asimov_paramset, llh_paramset=llh_paramset ) - ndim = len(mcmc_paramset) + ndim = len(llh_paramset) if args.mcmc_seed_type == MCMCSeedType.UNIFORM: p0 = mcmc_utils.flat_seed( - mcmc_paramset, nwalkers=args.nwalkers + llh_paramset, nwalkers=args.nwalkers ) elif args.mcmc_seed_type == MCMCSeedType.GAUSSIAN: p0 = mcmc_utils.gaussian_seed( - mcmc_paramset, nwalkers=args.nwalkers + llh_paramset, nwalkers=args.nwalkers ) samples = mcmc_utils.mcmc( @@ -320,275 +146,12 @@ def main(): mcmc_utils.save_chains(samples, outfile) plot_utils.chainer_plot( - infile = outfile+'.npy', - outfile = outfile[:5]+outfile[5:].replace('data', 'plots'), - outformat = ['pdf'], - args = args, - mcmc_paramset = mcmc_paramset - ) - - out = args.bayes_output+'/fr_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 - ) - if args.run_bayes_factor: - import pymultinest - - if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT: - fitter = gf_utils.setup_fitter(args, asimov_paramset) - else: fitter = None - - p = mcmc_paramset.from_tag(ParamTag.SCALE, invert=True) - n_params = len(p) - prior_ranges = p.seeds - - def CubePrior(cube, ndim, nparams): - # default are uniform priors - return ; - - arr = [] - for s_idx, sc in enumerate(scan_scales): - if args.bayes_eval_bin is not None: - if args.bayes_eval_bin == s_idx: - out += '_scale_{0:.0E}'.format(np.power(10, sc)) - else: continue - - print '== SCALE = {0:.0E}'.format(np.power(10, sc)) - theta = np.zeros(n_params) - def lnProb(cube, ndim, nparams): - for i in range(ndim): - prange = prior_ranges[i][1] - prior_ranges[i][0] - theta[i] = prange*cube[i] + prior_ranges[i][0] - theta_ = np.array(theta.tolist() + [sc]) - # print 'mcmc_paramset', mcmc_paramset - return llh_utils.triangle_llh( - theta=theta_, - args=args, - asimov_paramset=asimov_paramset, - mcmc_paramset=mcmc_paramset, - fitter=fitter - ) - - prefix = 'mnrun/' + os.path.basename(out) + '_' - if args.bayes_eval_bin is not None: - prefix += '{0}_'.format(s_idx) - print 'begin running evidence calculation for {0}'.format(prefix) - result = pymultinest.run( - LogLikelihood=lnProb, - Prior=CubePrior, - n_dims=n_params, - importance_nested_sampling=True, - n_live_points=args.bayes_live_points, - evidence_tolerance=args.bayes_tolerance, - outputfiles_basename=prefix, - resume=False, - verbose=True - ) - - analyzer = pymultinest.Analyzer( - outputfiles_basename=prefix, n_params=n_params - ) - a_lnZ = analyzer.get_stats()['global evidence'] - print 'Evidence = {0}'.format(a_lnZ) - arr.append([sc, a_lnZ]) - - misc_utils.make_dir(out) - np.save(out+'.npy', np.array(arr)) - - dirname = os.path.dirname(out) - plot_utils.bayes_factor_plot( - dirname=dirname, outfile=out, outformat=['png'], args=args - ) - - out = args.angles_lim_output+'/fr_an_evidence' + misc_utils.gen_identifier(args) - if args.run_angles_limit: - import pymultinest - - scenarios = [ - [np.sin(np.pi/2.)**2, 0, 0, 0], - [0, np.cos(np.pi/2.)**4, 0, 0], - [0, 0, np.sin(np.pi/2.)**2, 0], - ] - p = mcmc_paramset.from_tag([ParamTag.SCALE, ParamTag.MMANGLES], invert=True) - n_params = len(p) - prior_ranges = p.seeds - - if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT: - fitter = gf_utils.setup_fitter(args, asimov_paramset) - else: fitter = None - - def CubePrior(cube, ndim, nparams): - # default are uniform priors - return ; - - if args.bayes_eval_bin is not None: - data = np.zeros((len(scenarios), 1, 2)) - else: data = np.zeros((len(scenarios), args.bayes_bins, 2)) - mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES) - sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0] - for idx, scen in enumerate(scenarios): - scales, evidences = [], [] - for yidx, an in enumerate(mm_angles): - an.value = scen[yidx] - for s_idx, sc in enumerate(scan_scales): - if args.bayes_eval_bin is not None: - if args.bayes_eval_bin == s_idx: - if idx == 0: - out += '_scale_{0:.0E}'.format(np.power(10, sc)) - else: continue - - print '== SCALE = {0:.0E}'.format(np.power(10, sc)) - sc_angles.value = sc - def lnProb(cube, ndim, nparams): - for i in range(ndim): - prange = prior_ranges[i][1] - prior_ranges[i][0] - p[i].value = prange*cube[i] + prior_ranges[i][0] - for name in p.names: - mcmc_paramset[name].value = p[name].value - theta = mcmc_paramset.values - # print 'theta', theta - # print 'mcmc_paramset', mcmc_paramset - return llh_utils.triangle_llh( - theta=theta, - args=args, - asimov_paramset=asimov_paramset, - mcmc_paramset=mcmc_paramset, - fitter=fitter - ) - prefix = 'mnrun/' + os.path.basename(out) + '_' - if args.bayes_eval_bin is not None: - prefix += '{0}_{1}_'.format(idx, s_idx) - print 'begin running evidence calculation for {0}'.format(prefix) - result = pymultinest.run( - LogLikelihood=lnProb, - Prior=CubePrior, - n_dims=n_params, - importance_nested_sampling=True, - n_live_points=args.bayes_live_points, - evidence_tolerance=args.bayes_tolerance, - outputfiles_basename=prefix, - resume=False, - verbose=True - ) - - analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params) - a_lnZ = analyzer.get_stats()['global evidence'] - print 'Evidence = {0}'.format(a_lnZ) - scales.append(sc) - evidences.append(a_lnZ) - - for i, d in enumerate(evidences): - data[idx][i][0] = scales[i] - data[idx][i][1] = d - - misc_utils.make_dir(out) - print 'saving to {0}.npy'.format(out) - np.save(out+'.npy', np.array(data)) - - dirname = os.path.dirname(out) - plot_utils.plot_BSM_angles_limit( - dirname=dirname, outfile=outfile, outformat=['png'], - args=args, bayesian=True + infile = outfile+'.npy', + outfile = outfile[:5]+outfile[5:].replace('data', 'plots'), + outformat = ['pdf'], + args = args, + llh_paramset = llh_paramset ) - - out = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args) - if args.run_angles_correlation: - if args.bayes_eval_bin is None: assert 0 - import pymultinest - - scenarios = [ - [1, 0, 0, 0], - [0, 1, 0, 0], - [0, 0, 1, 0], - ] - nuisance = mcmc_paramset.from_tag(ParamTag.NUISANCE) - mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES) - sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0] - - if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT: - fitter = gf_utils.setup_fitter(args, asimov_paramset) - else: fitter = None - - def CubePrior(cube, ndim, nparams): - # default are uniform priors - return ; - - scan_angles = np.linspace(0, 1, args.bayes_bins) - - if args.bayes_eval_bin is not None: - data = np.zeros((len(scenarios), 1, 1, 3)) - else: data = np.zeros((len(scenarios), scale_bins, angle_bins, 3)) - for idx, scen in enumerate(scenarios): - for an in mm_angles: - an.value = 0 - keep = mcmc_paramset.from_tag(ParamTag.MMANGLES)[idx] - q = ParamSet(nuisance.params + [x for x in mm_angles if x.name != keep.name]) - n_params = len(q) - prior_ranges = q.seeds - - scales, angles, evidences = [], [], [] - for s_idx, sc in enumerate(scan_scales): - for a_idx, an in enumerate(scan_angles): - index = s_idx*args.bayes_bins + a_idx - if args.bayes_eval_bin is not None: - if args.bayes_eval_bin == index: - if idx == 0: - out += '_scale_{0:.0E}'.format(np.power(10, sc)) - else: continue - - print '== SCALE = {0:.0E}'.format(np.power(10, sc)) - print '== ANGLE = {0:.2f}'.format(an) - sc_angles.value = sc - keep.value = an - def lnProb(cube, ndim, nparams): - for i in range(ndim-1): - prange = prior_ranges[i][1] - prior_ranges[i][0] - q[i].value = prange*cube[i] + prior_ranges[i][0] - for name in q.names: - mcmc_paramset[name].value = q[name].value - theta = mcmc_paramset.values - # print 'theta', theta - # print 'mcmc_paramset', mcmc_paramset - return llh_utils.triangle_llh( - theta=theta, - args=args, - asimov_paramset=asimov_paramset, - mcmc_paramset=mcmc_paramset, - fitter=fitter - ) - prefix = 'mnrun/' + os.path.basename(out) + '_' - if args.bayes_eval_bin is not None: - prefix += '{0}_{1}_{2}'.format(idx, s_idx, a_idx) - - print 'begin running evidence calculation for {0}'.format(prefix) - result = pymultinest.run( - LogLikelihood=lnProb, - Prior=CubePrior, - n_dims=n_params, - importance_nested_sampling=True, - n_live_points=args.bayes_live_points, - evidence_tolerance=args.bayes_tolerance, - outputfiles_basename=prefix, - resume=False, - verbose=True - ) - - analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params) - a_lnZ = analyzer.get_stats()['global evidence'] - print 'Evidence = {0}'.format(a_lnZ) - scales.append(sc) - angles.append(an) - evidences.append(a_lnZ) - - for i, d in enumerate(evidences): - data[idx][i][i][0] = scales[i] - data[idx][i][i][1] = angles[i] - data[idx][i][i][2] = d - - misc_utils.make_dir(out) - print 'saving to {0}.npy'.format(out) - np.save(out+'.npy', np.array(data)) - print "DONE!" diff --git a/fr_edit.py b/fr_edit.py new file mode 100755 index 0000000..df3b43b --- /dev/null +++ b/fr_edit.py @@ -0,0 +1,581 @@ +#! /usr/bin/env python +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : March 17, 2018 + +""" +HESE BSM flavour ratio analysis script +""" + +from __future__ import absolute_import, division + +import os +import argparse +from functools import partial + +import numpy as np + +import GolemFitPy as gf + +from utils import gf as gf_utils +from utils import likelihood as llh_utils +from utils import mcmc as mcmc_utils +from utils import misc as misc_utils +from utils import plot as plot_utils +from utils.enums import EnergyDependance, Likelihood, MCMCSeedType, ParamTag +from utils.fr import MASS_EIGENVALUES, normalise_fr, fr_to_angles +from utils.misc import enum_parse, Param, ParamSet + + +def define_nuisance(): + """Define the nuisance parameters to scan over with default values, + ranges and sigma. + """ + tag = ParamTag.NUISANCE + nuisance = ParamSet( + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + ) + return nuisance + + +def nuisance_argparse(parser): + nuisance_paramset = define_nuisance() + for parm in nuisance_paramset: + parser.add_argument( + '--'+parm.name, type=float, default=parm.value, + help=parm.name+' to inject' + ) + + +def get_paramsets(args): + """Make the paramsets for generating the Asmimov MC sample and also running + the MCMC. + """ + asimov_paramset = [] + mcmc_paramset = [] + if args.likelihood == Likelihood.GOLEMFIT: + nuisance_paramset = define_nuisance() + asimov_paramset.extend(nuisance_paramset.params) + mcmc_paramset.extend(nuisance_paramset.params) + for parm in nuisance_paramset: + parm.value = args.__getattribute__(parm.name) + tag = ParamTag.BESTFIT + flavour_angles = fr_to_angles(args.measured_ratio) + asimov_paramset.extend([ + Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag), + Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag), + ]) + asimov_paramset = ParamSet(asimov_paramset) + + if not args.fix_mixing and not args.fix_mixing_almost: + tag = ParamTag.MMANGLES + e = 1e-10 + mcmc_paramset.extend([ + Param(name='s_12^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag), + Param(name='c_13^4', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag), + Param(name='s_23^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag), + Param(name='dcp', value=np.pi, ranges=[0., 2*np.pi], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag) + ]) + if args.fix_mixing_almost: + tag = ParamTag.MMANGLES + mcmc_paramset.extend([ + Param(name='s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag) + ]) + if not args.fix_scale: + tag = ParamTag.SCALE + mcmc_paramset.append( + Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, + tex=r'{\rm log}_{10}\Lambda'+plot_utils.get_units(args.dimension), tag=tag) + ) + if not args.fix_source_ratio: + tag = ParamTag.SRCANGLES + mcmc_paramset.extend([ + Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag), + Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag) + ]) + mcmc_paramset = ParamSet(mcmc_paramset) + return asimov_paramset, mcmc_paramset + + +def process_args(args): + """Process the input args.""" + if args.fix_mixing and args.fix_scale: + raise NotImplementedError('Fixed mixing and scale not implemented') + if args.fix_mixing and args.fix_mixing_almost: + raise NotImplementedError( + '--fix-mixing and --fix-mixing-almost cannot be used together' + ) + if args.run_bayes_factor and args.fix_scale: + raise NotImplementedError( + '--run-bayes-factor and --fix-scale cannot be used together' + ) + + args.measured_ratio = normalise_fr(args.measured_ratio) + if args.fix_source_ratio: + args.source_ratio = normalise_fr(args.source_ratio) + + if args.energy_dependance is EnergyDependance.SPECTRAL: + args.binning = np.logspace( + np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1 + ) + + if not args.fix_scale: + if args.energy_dependance is EnergyDependance.MONO: + args.scale = np.power( + 10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \ + np.log10(args.energy**(args.dimension-3)) + 6 + ) + elif args.energy_dependance is EnergyDependance.SPECTRAL: + args.scale = np.power( + 10, np.round( + np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \ + - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3)) + ) + ) + """Estimate the scale of when to expect the BSM term to contribute.""" + args.scale_region = (args.scale/args.scale_region, args.scale*args.scale_region) + + if args.bayes_eval_bin.lower() == 'all': + args.bayes_eval_bin = None + else: + args.bayes_eval_bin = map(int, args.bayes_eval_bin.split()) + + +def parse_args(args=None): + """Parse command line arguments""" + parser = argparse.ArgumentParser( + description="BSM flavour ratio analysis", + formatter_class=misc_utils.SortingHelpFormatter, + ) + parser.add_argument( + '--measured-ratio', type=int, nargs=3, default=[1, 1, 1], + help='Set the central value for the measured flavour ratio at IceCube' + ) + parser.add_argument( + '--sigma-ratio', type=float, default=0.01, + help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH' + ) + parser.add_argument( + '--fix-source-ratio', type=misc_utils.parse_bool, default='False', + help='Fix the source flavour ratio' + ) + parser.add_argument( + '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance), + choices=EnergyDependance, + help='Type of energy dependance to use in the BSM fit' + ) + parser.add_argument( + '--spectral-index', default=-2, type=int, + help='Spectral index for spectral energy dependance' + ) + parser.add_argument( + '--binning', default=[1e4, 1e7, 5], type=float, nargs=3, + help='Binning for spectral energy dependance' + ) + parser.add_argument( + '--run-bayes-factor', type=misc_utils.parse_bool, default='False', + help='Make the bayes factor plot for the scale' + ) + parser.add_argument( + '--run-angles-limit', type=misc_utils.parse_bool, default='False', + help='Make the limit vs BSM angles plot' + ) + parser.add_argument( + '--run-angles-correlation', type=misc_utils.parse_bool, default='False', + help='Make the limit vs BSM angles plot' + ) + parser.add_argument( + '--bayes-bins', type=int, default=10, + help='Number of bins for the Bayes factor plot' + ) + parser.add_argument( + '--bayes-eval-bin', type=str, default='all', + help='Which bin to evalaute for Bayes factor plot' + ) + parser.add_argument( + '--bayes-live-points', type=int, default=400, + help='Number of live points for MultiNest evaluations' + ) + parser.add_argument( + '--bayes-tolerance', type=float, default=0.5, + help='Tolerance for MultiNest' + ) + parser.add_argument( + '--bayes-output', type=str, default='./mnrun/', + help='Folder to store MultiNest evaluations' + ) + parser.add_argument( + '--angles-lim-output', type=str, default='./mnrun/', + help='Folder to store MultiNest evaluations' + ) + parser.add_argument( + '--angles-corr-output', type=str, default='./mnrun/', + help='Folder to store MultiNest evaluations' + ) + parser.add_argument( + '--source-ratio', type=int, nargs=3, default=[2, 1, 0], + help='Set the source flavour ratio for the case when you want to fix it' + ) + parser.add_argument( + '--no-bsm', type=misc_utils.parse_bool, default='False', + help='Turn off BSM terms' + ) + parser.add_argument( + '--fix-mixing', type=misc_utils.parse_bool, default='False', + help='Fix all mixing parameters to bi-maximal mixing' + ) + parser.add_argument( + '--fix-mixing-almost', type=misc_utils.parse_bool, default='False', + help='Fix all mixing parameters except s23' + ) + parser.add_argument( + '--fix-scale', type=misc_utils.parse_bool, default='False', + help='Fix the new physics scale' + ) + parser.add_argument( + '--scale', type=float, default=1e-30, + help='Set the new physics scale' + ) + parser.add_argument( + '--scale-region', type=float, default=1e10, + help='Set the size of the box to scan for the scale' + ) + parser.add_argument( + '--dimension', type=int, default=3, + help='Set the new physics dimension to consider' + ) + parser.add_argument( + '--energy', type=float, default=1000, + help='Set the energy scale' + ) + parser.add_argument( + '--seed', type=int, default=99, + help='Set the random seed value' + ) + parser.add_argument( + '--threads', type=misc_utils.thread_type, default='1', + help='Set the number of threads to use (int or "max")' + ) + parser.add_argument( + '--outfile', type=str, default='./untitled', + help='Path to output chains' + ) + llh_utils.likelihood_argparse(parser) + mcmc_utils.mcmc_argparse(parser) + gf_utils.gf_argparse(parser) + plot_utils.plot_argparse(parser) + nuisance_argparse(parser) + if args is None: return parser.parse_args() + else: return parser.parse_args(args.split()) + + +def main(): + args = parse_args() + process_args(args) + misc_utils.print_args(args) + + np.random.seed(args.seed) + + asimov_paramset, mcmc_paramset = get_paramsets(args) + outfile = misc_utils.gen_outfile_name(args) + print '== {0:<25} = {1}'.format('outfile', outfile) + + if args.run_mcmc: + if args.likelihood is Likelihood.GOLEMFIT: + fitter = gf_utils.setup_fitter(args, asimov_paramset) + else: + fitter = None + + ln_prob = partial( + llh_utils.ln_prob, args=args, fitter=fitter, + asimov_paramset=asimov_paramset, mcmc_paramset=mcmc_paramset + ) + + ndim = len(mcmc_paramset) + if args.mcmc_seed_type == MCMCSeedType.UNIFORM: + p0 = mcmc_utils.flat_seed( + mcmc_paramset, nwalkers=args.nwalkers + ) + elif args.mcmc_seed_type == MCMCSeedType.GAUSSIAN: + p0 = mcmc_utils.gaussian_seed( + mcmc_paramset, nwalkers=args.nwalkers + ) + + samples = mcmc_utils.mcmc( + p0 = p0, + ln_prob = ln_prob, + ndim = ndim, + nwalkers = args.nwalkers, + burnin = args.burnin, + nsteps = args.nsteps, + threads = 1 + # TODO(shivesh): broken because you cannot pickle a GolemFitPy object + # threads = misc_utils.thread_factors(args.threads)[0] + ) + mcmc_utils.save_chains(samples, outfile) + + plot_utils.chainer_plot( + infile = outfile+'.npy', + outfile = outfile[:5]+outfile[5:].replace('data', 'plots'), + outformat = ['pdf'], + args = args, + mcmc_paramset = mcmc_paramset + ) + + out = args.bayes_output+'/fr_fr_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 + ) + if args.run_bayes_factor: + from scipy.optimize import minimize + def fit_flags(flag_dict): + flags = gf.FitParametersFlag() + for key in flag_dict.iterkeys(): + flags.__setattr__(key, flag_dict[key]) + return flags + + args.likelihood = Likelihood.GF_FREQ + default_flags = { + # False means it's not fixed in minimization + 'astroFlavorAngle1' : True, + 'astroFlavorAngle2' : True, + 'astroENorm' : False, + 'astroMuNorm' : False, + 'astroTauNorm' : False, + 'convNorm' : False, + 'promptNorm' : False, + 'muonNorm' : False, + 'astroNorm' : False, + 'astroParticleBalance' : True, + 'astroDeltaGamma' : False, + 'cutoffEnergy' : True, + 'CRDeltaGamma' : True, + 'piKRatio' : True, + 'NeutrinoAntineutrinoRatio' : True, + 'darkNorm' : True, + 'domEfficiency' : True, + 'holeiceForward' : True, + 'anisotropyScale' : True, + 'astroNormSec' : True, + 'astroDeltaGammaSec' : True + } + + mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True) + scale = mcmc_paramset.from_tag(ParamTag.SCALE) + + if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ: + fitter = gf_utils.setup_fitter(args, asimov_paramset) + fitter.SetFitParametersFlag(fit_flags(default_flags)) + else: fitter = None + + data = np.zeros((args.bayes_bins, 2)) + for s_idx, sc in enumerate(scan_scales): + print '== SCALE = {0:.0E}'.format(np.power(10, sc)) + def fn(scen): + theta = map(float, scen) + [0, sc] + try: + llh = llh_utils.triangle_llh( + theta=theta, args=args, asimov_paramset=asimov_paramset, + mcmc_paramset=mcmc_paramset_freq, fitter=fitter + ) + print 'llh', llh + except: + print 'forbidden' + return 9e99 + return -llh + res = minimize(fun=fn, x0=[0.1, 0.1, 0.1], method='L-BFGS-B', + bounds=[(0, 1), (0, 1), (0, 1)]) + llh = fn(res.x) + data[s_idx][0] = sc + data[s_idx][1] = llh + + misc_utils.make_dir(out) + np.save(out+'.npy', np.array(data)) + + # dirname = os.path.dirname(out) + # plot_utils.bayes_factor_plot( + # dirname=dirname, outfile=out, outformat=['png'], args=args + # ) + + out = args.angles_lim_output+'/fr_anfr_evidence' + misc_utils.gen_identifier(args) + if args.run_angles_limit: + def fit_flags(flag_dict): + flags = gf.FitParametersFlag() + for key in flag_dict.iterkeys(): + flags.__setattr__(key, flag_dict[key]) + return flags + + args.likelihood = Likelihood.GF_FREQ + default_flags = { + # False means it's not fixed in minimization + 'astroFlavorAngle1' : True, + 'astroFlavorAngle2' : True, + # 'astroENorm' : True, + # 'astroMuNorm' : True, + # 'astroTauNorm' : True, + 'convNorm' : False, + 'promptNorm' : False, + 'muonNorm' : False, + 'astroNorm' : False, + 'astroParticleBalance' : True, + 'astroDeltaGamma' : False, + 'cutoffEnergy' : True, + 'CRDeltaGamma' : True, + 'piKRatio' : True, + 'NeutrinoAntineutrinoRatio' : True, + 'darkNorm' : True, + 'domEfficiency' : True, + 'holeiceForward' : True, + 'anisotropyScale' : True, + 'astroNormSec' : True, + 'astroDeltaGammaSec' : True + } + + mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True) + scenarios = [ + [np.sin(np.pi/2.)**2, 0, 0, 0], + [0, np.cos(np.pi/2.)**4, 0, 0], + [0, 0, np.sin(np.pi/2.)**2, 0], + ] + + if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ: + fitter = gf_utils.setup_fitter(args, asimov_paramset) + fitter.SetFitParametersFlag(fit_flags(default_flags)) + else: fitter = None + + data = np.zeros((len(scenarios), args.bayes_bins, 2)) + mm_angles = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES) + sc_angles = mcmc_paramset_freq.from_tag(ParamTag.SCALE)[0] + for idx, scen in enumerate(scenarios): + scales, llhs = [], [] + for yidx, an in enumerate(mm_angles): + an.value = scen[yidx] + for s_idx, sc in enumerate(scan_scales): + print '== SCALE = {0:.0E}'.format(np.power(10, sc)) + theta = scen + [sc] + print 'theta', theta + llh = llh_utils.triangle_llh( + theta=theta, args=args, asimov_paramset=asimov_paramset, + mcmc_paramset=mcmc_paramset_freq, fitter=fitter + ) + print 'llh', llh + scales.append(sc) + llhs.append(llh) + + for i, d in enumerate(llhs): + data[idx][i][0] = scales[i] + data[idx][i][1] = d + + misc_utils.make_dir(out) + print 'saving to {0}.npy'.format(out) + np.save(out+'.npy', np.array(data)) + + dirname = os.path.dirname(out) + plot_utils.plot_BSM_angles_limit( + dirname=dirname, outfile=outfile, outformat=['png'], + args=args, bayesian=True + ) + + out = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args) + if args.run_angles_correlation: + def fit_flags(flag_dict): + flags = gf.FitParametersFlag() + for key in flag_dict.iterkeys(): + flags.__setattr__(key, flag_dict[key]) + return flags + + args.likelihood = Likelihood.GF_FREQ + default_flags = { + # False means it's not fixed in minimization + 'astroFlavorAngle1' : True, + 'astroFlavorAngle2' : True, + # 'astroENorm' : True, + # 'astroMuNorm' : True, + # 'astroTauNorm' : True, + 'convNorm' : False, + 'promptNorm' : False, + 'muonNorm' : False, + 'astroNorm' : False, + 'astroParticleBalance' : True, + 'astroDeltaGamma' : False, + 'cutoffEnergy' : True, + 'CRDeltaGamma' : True, + 'piKRatio' : True, + 'NeutrinoAntineutrinoRatio' : True, + 'darkNorm' : True, + 'domEfficiency' : True, + 'holeiceForward' : True, + 'anisotropyScale' : True, + 'astroNormSec' : True, + 'astroDeltaGammaSec' : True + } + + mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True) + + scenarios = [ + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + ] + mm_angles = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES) + sc_angles = mcmc_paramset_freq.from_tag(ParamTag.SCALE)[0] + + if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ: + fitter = gf_utils.setup_fitter(args, asimov_paramset) + fitter.SetFitParametersFlag(fit_flags(default_flags)) + else: fitter = None + + scan_angles = np.linspace(0, 1, args.bayes_bins) + + data = np.zeros((len(scenarios), args.bayes_bins, args.bayes_bins, 3)) + for idx, scen in enumerate(scenarios): + for an in mm_angles: + an.value = 0 + keep = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES)[idx] + + for s_idx, sc in enumerate(scan_scales): + scales, angles, llhs = [], [], [] + for a_idx, an in enumerate(scan_angles): + print '== SCALE = {0:.0E}'.format(np.power(10, sc)) + print '== ANGLE = {0:.2f}'.format(an) + sc_angles.value = sc + s2_2 = an + if s_idx == 0: keep.value = np.sin(np.arcsin(np.sqrt(s2_2))/2.)**2 + elif s_idx == 1: keep.value = np.cos(np.arcsin(np.sqrt(s2_2))/2.)**4 + elif s_idx == 2: keep.value = np.sin(np.arcsin(np.sqrt(s2_2))/2.)**2 + theta = mcmc_paramset_freq.values + print 'mcmc_paramset_freq', mcmc_paramset_freq + try: + llh = llh_utils.triangle_llh( + theta=theta, args=args, asimov_paramset=asimov_paramset, + mcmc_paramset=mcmc_paramset_freq, fitter=fitter + ) + except AssertionError: + print 'forbidden by unitarity' + data[idx][s_idx][a_idx][0] = np.nan + data[idx][s_idx][a_idx][1] = np.nan + data[idx][s_idx][a_idx][2] = np.nan + continue + print 'llh', llh + data[idx][s_idx][a_idx][0] = sc + data[idx][s_idx][a_idx][1] = an + data[idx][s_idx][a_idx][2] = llh + + misc_utils.make_dir(out) + print 'saving to {0}.npy'.format(out) + np.save(out+'.npy', np.array(data)) + + print "DONE!" + + +main.__doc__ = __doc__ + + +if __name__ == '__main__': + main() + diff --git a/plot_bf.py b/plot_bf.py index faefcde..13664b9 100644..100755 --- a/plot_bf.py +++ b/plot_bf.py @@ -24,7 +24,7 @@ from matplotlib.offsetbox import AnchoredText import fr from utils import misc as misc_utils from utils.fr import normalise_fr -from utils.plot import myround, get_units +from utils.plot import bayes_factor_plot, myround, get_units rc('text', usetex=False) @@ -33,11 +33,12 @@ 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), + (1, 1, 1, 0, 1, 0), ] # FR -dimension = [3, 4, 5, 6, 7, 8] +dimension = [3, 6] +# dimension = [3, 4, 5, 6, 7, 8] sigma_ratio = ['0.01'] energy_dependance = 'spectral' spectral_index = -2 @@ -68,15 +69,16 @@ bayes_eval_bin = True # set to 'all' to run normally # Plot plot_bayes = False -plot_angles_limit = False -plot_angles_corr = True +plot_angles_limit = True +plot_angles_corr = False outformat = ['png'] -significance = np.log(10**(1/2.)) -# significance = np.log(10**(3/2.)) +# 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) @@ -96,7 +98,7 @@ for i_dim, dim in enumerate(dimension): 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'.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) + 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) @@ -110,18 +112,33 @@ for i_dim, dim in enumerate(dimension): 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 - for sc in scan_scales: + 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,:] - if plot_angles_corr: - # TODO(shivesh) - assert 0 - if len(lf.shape) == 3: lf = lf[:,0,:] raw.append(lf) except IOError: fail += 1 @@ -141,13 +158,28 @@ for i_dim, dim in enumerate(dimension): 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: - # TODO(shivesh) 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)) @@ -160,7 +192,7 @@ if plot_bayes: 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 / GeV^{-d+4}$') + 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 @@ -187,8 +219,11 @@ if plot_bayes: else: print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, null) - yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) - ax.set_ylim(yranges) + # 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(): @@ -197,7 +232,7 @@ if plot_bayes: ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) for of in outformat: - fig.savefig('./images/bayes_factor.'+of, bbox_inches='tight', dpi=150) + 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'} @@ -210,15 +245,16 @@ if plot_angles_limit: ax.set_xlim(0, len(xticks)+1) ax.set_xticklabels([''] + xticks + ['']) ax.set_xlabel(r'BSM operator angle') - ylabel = r'${\rm log}_{10} \Lambda' + get_units(dim) + r'$' + 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[0] + 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] @@ -249,4 +285,130 @@ if plot_angles_limit: ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1) for of in outformat: - fig.savefig('./images/angles_limit_DIM{0}'.format(dim)+'.'+of, bbox_inches='tight', dpi=150) + 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) + @@ -0,0 +1,480 @@ +#! /usr/bin/env python +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : March 17, 2018 + +""" +HESE BSM flavour ratio analysis script +""" + +from __future__ import absolute_import, division + +import os +import argparse +from functools import partial + +import numpy as np +import numpy.ma as ma + +from utils import fr as fr_utils +from utils import gf as gf_utils +from utils import likelihood as llh_utils +from utils import misc as misc_utils +from utils import plot as plot_utils +from utils.enums import EnergyDependance, Likelihood, ParamTag +from utils.enums import SensitivityCateg, StatCateg +from utils.fr import estimate_scale, normalise_fr +from utils.misc import enum_parse, parse_bool +from utils.param import Param, ParamSet, get_paramsets + +from utils import multinest as mn_utils + + +def define_nuisance(): + """Define the nuisance parameters to scan over with default values, + ranges and sigma. + """ + tag = ParamTag.NUISANCE + nuisance = ParamSet( + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + ) + return nuisance + + +def nuisance_argparse(parser): + nuisance_paramset = define_nuisance() + for parm in nuisance_paramset: + parser.add_argument( + '--'+parm.name, type=float, default=parm.value, + help=parm.name+' to inject' + ) + + +def process_args(args): + """Process the input args.""" + if args.fix_mixing and args.fix_scale: + raise NotImplementedError('Fixed mixing and scale not implemented') + if args.fix_mixing and args.fix_mixing_almost: + raise NotImplementedError( + '--fix-mixing and --fix-mixing-almost cannot be used together' + ) + if args.mn_run and args.fix_scale: + raise NotImplementedError( + '--mn-run and --fix-scale cannot be used together' + ) + + args.measured_ratio = normalise_fr(args.measured_ratio) + if args.fix_source_ratio: + args.source_ratio = normalise_fr(args.source_ratio) + + if args.energy_dependance is EnergyDependance.SPECTRAL: + args.binning = np.logspace( + np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1 + ) + + if not args.fix_scale: + args.scale = estimate_scale(args) + args.scale_region = (args.scale/args.scale_region, args.scale*args.scale_region) + + if args.mn_eval_bin.lower() == 'all': + args.mn_eval_bin = None + else: + args.mn_eval_bin = int(args.mn_eval_bin) + + +def parse_args(args=None): + """Parse command line arguments""" + parser = argparse.ArgumentParser( + description="BSM flavour ratio analysis", + formatter_class=misc_utils.SortingHelpFormatter, + ) + parser.add_argument( + '--seed', type=misc_utils.seed_parse, default='25', + help='Set the random seed value' + ) + parser.add_argument( + '--threads', type=misc_utils.thread_type, default='1', + help='Set the number of threads to use (int or "max")' + ) + parser.add_argument( + '--outfile', type=str, default='./untitled', + help='Path to output chains' + ) + parser.add_argument( + '--run-method', default='full', + type=partial(enum_parse, c=SensitivityCateg), choices=SensitivityCateg, + help='Choose which type of sensivity plot to make' + ) + parser.add_argument( + '--stat-method', default='bayesian', + type=partial(enum_parse, c=StatCateg), choices=StatCateg, + help='Statistical method to employ' + ) + fr_utils.fr_argparse(parser) + gf_utils.gf_argparse(parser) + llh_utils.likelihood_argparse(parser) + mn_utils.mn_argparse(parser) + nuisance_argparse(parser) + if args is None: return parser.parse_args() + else: return parser.parse_args(args.split()) + +def main(): + args = parse_args() + process_args(args) + misc_utils.print_args(args) + + if args.seed is not None: + np.random.seed(args.seed) + + asimov_paramset, llh_paramset = get_paramsets(args) + scale = llh_paramset.from_tag(ParamTag.SCALE)[0] + outfile = misc_utils.gen_outfile_name(args) + print '== {0:<25} = {1}'.format('outfile', outfile) + + mmangles = llh_paramset.from_tag(ParamTag.MMANGLES) + if args.run_method is SensitivityCateg.FULL: + mn_paramset_arr = [llh_paramset.from_tag(ParamTag.SCALE, invert=True)] + elif args.run_method is SensitivityCateg.FIXED_ANGLE or \ + args.run_method is SensitivityCateg.CORR_ANGLE: + nscale_pmset = llh_paramset.from_tag(ParamTag.SCALE, invert=True) + mn_paramset_arr = [] + for x in xrange(3): + mn_paramset_arr.append( + ParamSet([prms for prms in nscale_pmset + if mmangles[x].name != prms.name]) + ) + + out = args.outfile+'/{0}/{1}/fr_evidence'.format(args.stat_method, args.run_method) \ + + misc_utils.gen_identifier(args) + if args.mn_run: + if args.likelihood is Likelihood.GOLEMFIT: + fitter = gf_utils.setup_fitter(args, asimov_paramset) + else: fitter = None + + scan_scales = np.linspace( + np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.mn_bins + ) + if args.run_method is SensitivityCateg.CORR_ANGLE: + scan_angles = np.linspace(0, 1, eval_dim) + else: scan_angles = np.array([0]) + print 'scan_scales', scan_scales + print 'scan_angles', scan_angles + + if args.mn_eval_bin is None: + eval_dim = args.mn_bins + else: eval_dim = 1 + + if args.run_method is SensitivityCateg.FULL: + evidence_arr = np.full((eval_dim, 2), np.nan) + elif args.run_method is SensitivityCateg.FIXED_ANGLE: + evidence_arr = np.full((len(mn_paramset_arr), eval_dim, 2), np.nan) + elif args.run_method is SensitivityCateg.CORR_ANGLE: + evidence_arr = np.full((len(mn_paramset_arr), eval_dim, eval_dim, 3), np.nan) + + + for idx_scen, mn_paramset in enumerate(mn_paramset_arr): + print '|||| SCENARIO = {0}'.format(idx_scen) + for x in mmangles: x.value = 0. + if args.run_method is SensitivityCateg.FIXED_ANGLE: + if idx_scen == 0 or idx_scen == 2: + mmangles[idx_scen].value = np.sin(np.pi/2.)**2 + """s_12^2 or s_23^2""" + elif idx_scen == 1: + mmangles[idx_scen].value = np.cos(np.pi/2.)**4 + """c_13^4""" + + for idx_an, an in enumerate(scan_angles): + if args.run_method is SensitivityCateg.CORR_ANGLE: + print '|||| ANGLE = {0:<04.2}'.format(an) + mmangles[idx_an].value = an + + for idx_sc, sc in enumerate(scan_scales): + if args.mn_eval_bin is not None: + if idx_sc == args.mn_eval_bin: + out += '_scale_{0:.0E}'.format(np.power(10, sc)) + if args.run_method is SensitivityCateg.CORR_ANGLE: + out += '_angle_{0:>03}'.format(int(an*100)) + break + else: continue + + print '|||| SCALE = {0:.0E}'.format(np.power(10, sc)) + scale.value = sc + try: + a_lnZ = mn_utils.mn_evidence( + mn_paramset = mn_paramset, + llh_paramset = llh_paramset, + asimov_paramset = asimov_paramset, + args = args, + fitter = fitter + ) + except: + print 'Failed run, continuing' + continue + print '## Evidence = {0}'.format(a_lnZ) + if args.run_method is SensitivityCateg.FULL: + evidence_arr[idx_sc] = np.array([sc, a_lnZ]) + elif args.run_method is SensitivityCateg.FIXED_ANGLE: + evidence_arr[idx_scen][idx_sc] = np.array([sc, a_lnZ]) + elif args.run_method is SensitivityCateg.CORR_ANGLE: + evidence_arr[idx_scen][idx_an][idx_sc] = np.array([an, sc, a_lnZ]) + + misc_utils.make_dir(out) + print 'Saving to {0}'.format(out+'.npy') + np.save(out+'.npy', np.array(evidence_arr)) + + if args.plot_multinest: + if args.mn_run: raw = evidence_arr + else: raw = np.load(out+'.npy') + data = ma.masked_invalid(raw, 0) + + basename = os.path.dirname(out) + '/mnrun/' + os.path.basename(out) + baseoutfile = basename[:5]+basename[5:].replace('data', 'plots') + if args.run_method is SensitivityCateg.FULL: + plot_utils.plot_multinest( + data = data, + outfile = baseoutfile, + outformat = ['png'], + args = args, + scale_param = scale + ) + elif args.run_method is SensitivityCateg.FIXED_ANGLE: + for idx_scen in xrange(len(mn_paramset_arr)): + print '|||| SCENARIO = {0}'.format(idx_scen) + outfile = baseoutfile + '_SCEN{0}'.format(idx_scen) + if idx_scen == 0: label = r'$\mathcal{O}_{12}=\frac{1}{2}$' + elif idx_scen == 1: label = r'$\mathcal{O}_{13}=\frac{1}{2}$' + elif idx_scen == 2: label = r'$\mathcal{O}_{23}=\frac{1}{2}$' + plot_utils.plot_multinest( + data = data[idx_scen], + outfile = outfile, + outformat = ['png'], + args = args, + scale_param = scale, + label = label + ) + elif args.run_method is SensitivityCateg.CORR_ANGLE: + for idx_scen in xrange(len(mn_paramset_arr)): + print '|||| SCENARIO = {0}'.format(idx_scen) + basescenoutfile = baseoutfile + '_SCEN{0}'.format(idx_scen) + if idx_scen == 0: label = r'$\mathcal{O}_{12}=' + elif idx_scen == 1: label = r'$\mathcal{O}_{13}=' + elif idx_scen == 2: label = r'$\mathcal{O}_{23}=' + for idx_an, an in enumerate(scan_angles): + print '|||| ANGLE = {0:<04.2}'.format(an) + outfile = basescenoutfile + '_ANGLE{0}'.format(idx_an) + label += r'{0:<04.2}$'.format(an) + plot_utils.plot_multinest( + data = data[idx_scen][idx_an][:,1:], + outfile = outfile, + outformat = ['png'], + args = args, + scale_param = scale, + label = label + ) + + # dirname = os.path.dirname(out) + # plot_utils.bayes_factor_plot( + # dirname=dirname, outfile=out, outformat=['png'], args=args + # ) + + # out = args.angles_lim_output+'/fr_an_evidence' + misc_utils.gen_identifier(args) + # if args.run_angles_limit: + # import pymultinest + + # scenarios = [ + # [np.sin(np.pi/2.)**2, 0, 0, 0], + # [0, np.cos(np.pi/2.)**4, 0, 0], + # [0, 0, np.sin(np.pi/2.)**2, 0], + # ] + # p = llh_paramset.from_tag([ParamTag.SCALE, ParamTag.MMANGLES], invert=True) + # n_params = len(p) + # prior_ranges = p.seeds + + # if not args.run_llh and args.likelihood is Likelihood.GOLEMFIT: + # fitter = gf_utils.setup_fitter(args, asimov_paramset) + # else: fitter = None + + # def CubePrior(cube, ndim, nparams): + # # default are uniform priors + # return ; + + # if args.bayes_eval_bin is not None: + # data = np.zeros((len(scenarios), 1, 2)) + # else: data = np.zeros((len(scenarios), args.bayes_bins, 2)) + # mm_angles = llh_paramset.from_tag(ParamTag.MMANGLES) + # sc_angles = llh_paramset.from_tag(ParamTag.SCALE)[0] + # for idx, scen in enumerate(scenarios): + # scales, evidences = [], [] + # for yidx, an in enumerate(mm_angles): + # an.value = scen[yidx] + # for s_idx, sc in enumerate(scan_scales): + # if args.bayes_eval_bin is not None: + # if s_idx in args.bayes_eval_bin: + # if idx == 0: + # out += '_scale_{0:.0E}'.format(np.power(10, sc)) + # else: continue + + # print '== SCALE = {0:.0E}'.format(np.power(10, sc)) + # sc_angles.value = sc + # def lnProb(cube, ndim, nparams): + # for i in range(ndim): + # prange = prior_ranges[i][1] - prior_ranges[i][0] + # p[i].value = prange*cube[i] + prior_ranges[i][0] + # for name in p.names: + # mcmc_paramset[name].value = p[name].value + # theta = mcmc_paramset.values + # # print 'theta', theta + # # print 'mcmc_paramset', mcmc_paramset + # return llh_utils.triangle_llh( + # theta=theta, + # args=args, + # asimov_paramset=asimov_paramset, + # mcmc_paramset=mcmc_paramset, + # fitter=fitter + # ) + # prefix = 'mnrun/' + os.path.basename(out) + '_' + # if args.bayes_eval_bin is not None: + # prefix += '{0}_{1}_'.format(idx, s_idx) + # print 'begin running evidence calculation for {0}'.format(prefix) + # result = pymultinest.run( + # LogLikelihood=lnProb, + # Prior=CubePrior, + # n_dims=n_params, + # importance_nested_sampling=True, + # n_live_points=args.bayes_live_points, + # evidence_tolerance=args.bayes_tolerance, + # outputfiles_basename=prefix, + # resume=False, + # verbose=True + # ) + + # analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params) + # a_lnZ = analyzer.get_stats()['global evidence'] + # print 'Evidence = {0}'.format(a_lnZ) + # scales.append(sc) + # evidences.append(a_lnZ) + + # for i, d in enumerate(evidences): + # data[idx][i][0] = scales[i] + # data[idx][i][1] = d + + # misc_utils.make_dir(out) + # print 'saving to {0}.npy'.format(out) + # np.save(out+'.npy', np.array(data)) + + # dirname = os.path.dirname(out) + # plot_utils.plot_BSM_angles_limit( + # dirname=dirname, outfile=outfile, outformat=['png'], + # args=args, bayesian=True + # ) + + # out = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args) + # if args.run_angles_correlation: + # if args.bayes_eval_bin is None: assert 0 + # import pymultinest + + # scenarios = [ + # [1, 0, 0, 0], + # [0, 1, 0, 0], + # [0, 0, 1, 0], + # ] + # nuisance = mcmc_paramset.from_tag(ParamTag.NUISANCE) + # mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES) + # sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0] + + # if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT: + # fitter = gf_utils.setup_fitter(args, asimov_paramset) + # else: fitter = None + + # def CubePrior(cube, ndim, nparams): + # # default are uniform priors + # return ; + + # scan_angles = np.linspace(0, 1, args.bayes_bins) + + # if args.bayes_eval_bin is not None: + # data = np.zeros((len(scenarios), 1, 1, 3)) + # else: data = np.zeros((len(scenarios), args.bayes_bins, args.bayes_bins, 3)) + # for idx, scen in enumerate(scenarios): + # for an in mm_angles: + # an.value = 0 + # keep = mcmc_paramset.from_tag(ParamTag.MMANGLES)[idx] + # q = ParamSet(nuisance.params + [x for x in mm_angles if x.name != keep.name]) + # n_params = len(q) + # prior_ranges = q.seeds + + # scales, angles, evidences = [], [], [] + # for s_idx, sc in enumerate(scan_scales): + # for a_idx, an in enumerate(scan_angles): + # index = s_idx*args.bayes_bins + a_idx + # if args.bayes_eval_bin is not None: + # if index in args.bayes_eval_bin: + # if idx == 0: + # out += '_idx_{0}'.format(index) + # else: continue + + # print '== SCALE = {0:.0E}'.format(np.power(10, sc)) + # print '== ANGLE = {0:.2f}'.format(an) + # sc_angles.value = sc + # keep.value = an + # def lnProb(cube, ndim, nparams): + # for i in range(ndim-1): + # prange = prior_ranges[i][1] - prior_ranges[i][0] + # q[i].value = prange*cube[i] + prior_ranges[i][0] + # for name in q.names: + # mcmc_paramset[name].value = q[name].value + # theta = mcmc_paramset.values + # # print 'theta', theta + # # print 'mcmc_paramset', mcmc_paramset + # return llh_utils.triangle_llh( + # theta=theta, + # args=args, + # asimov_paramset=asimov_paramset, + # mcmc_paramset=mcmc_paramset, + # fitter=fitter + # ) + # prefix = 'mnrun/' + os.path.basename(out) + '_' + # if args.bayes_eval_bin is not None: + # prefix += '{0}_{1}_{2}'.format(idx, s_idx, a_idx) + + # print 'begin running evidence calculation for {0}'.format(prefix) + # result = pymultinest.run( + # LogLikelihood=lnProb, + # Prior=CubePrior, + # n_dims=n_params, + # importance_nested_sampling=True, + # n_live_points=args.bayes_live_points, + # evidence_tolerance=args.bayes_tolerance, + # outputfiles_basename=prefix, + # resume=False, + # verbose=True + # ) + + # analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params) + # a_lnZ = analyzer.get_stats()['global evidence'] + # print 'Evidence = {0}'.format(a_lnZ) + # scales.append(sc) + # angles.append(an) + # evidences.append(a_lnZ) + + # for i, d in enumerate(evidences): + # data[idx][i][i][0] = scales[i] + # data[idx][i][i][1] = angles[i] + # data[idx][i][i][2] = d + + # misc_utils.make_dir(out) + # print 'saving to {0}.npy'.format(out) + # np.save(out+'.npy', np.array(data)) + + +main.__doc__ = __doc__ + + +if __name__ == '__main__': + main() + diff --git a/submitter/make_dag.py b/submitter/make_dag.py index 78b7bff..c6e7400 100644 --- a/submitter/make_dag.py +++ b/submitter/make_dag.py @@ -19,7 +19,7 @@ fix_sfr_mfr = [ (1, 1, 1, 1, 2, 0), # (1, 1, 0, 1, 2, 0), # (1, 2, 0, 1, 2, 0), - (1, 1, 1, 1, 0, 0), + # (1, 1, 1, 1, 0, 0), # (1, 1, 0, 1, 0, 0), # (1, 0, 0, 1, 0, 0), (1, 1, 1, 0, 1, 0), @@ -34,12 +34,12 @@ run_mcmc = 'False' burnin = 2500 nsteps = 10000 nwalkers = 60 -seed = 24 +seed = 'None' threads = 4 mcmc_seed_type = 'uniform' # FR -dimension = [6] +dimension = [4, 5, 7, 8] energy = [1e6] no_bsm = 'False' sigma_ratio = ['0.01'] @@ -66,13 +66,13 @@ ast = 'p2_0' data = 'real' # Bayes Factor -run_bayes_factor = 'False' +run_bayes_factor = 'True' run_angles_limit = 'False' -run_angles_correlation = 'True' -bayes_bins = 20 -bayes_live_points = 1000 +run_angles_correlation = 'False' +bayes_bins = 100 +bayes_live_points = 3000 bayes_tolerance = 0.01 -bayes_eval_bin = 'None' # set to 'all' to run normally +bayes_eval_bin = 'all' # set to 'all' to run normally # Plot plot_angles = 'False' @@ -80,8 +80,8 @@ plot_elements = 'False' plot_bayes = 'False' plot_angles_limit = 'False' -# outfile = 'dagman_FR.submit'.format(dimension[0]) -outfile = 'dagman_FR_angles_correlation_DIM{0}.submit'.format(dimension[0]) +outfile = 'dagman_FR_freq_fullscan_otherdims.submit' +# outfile = 'dagman_FR_bayes_freq.submit' golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit' condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/submit.sub' @@ -106,6 +106,7 @@ with open(outfile, 'w') as f: bayes_output = 'None' angles_lim_output = 'None' + angles_corr_output = 'None' for sig in sigma_ratio: print 'sigma', sig for frs in fix_sfr_mfr: diff --git a/test/LV.out b/test/LV.out deleted file mode 100644 index e9f2b7e..0000000 --- a/test/LV.out +++ /dev/null @@ -1,74 +0,0 @@ -test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 2u, std::allocator<double> > already registered; second conversion method ignored. - import GolemFitPy as gf -test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored. - import GolemFitPy as gf -reset_steering: 1 -reset_data: 1 -reset_npp: 1 -GolemFit constructor: checking paths -Constructing DiffuseWeightMaker -Loading data -Loading HESE data. -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/data/HESEData.txt -Loaded HESE 80 events. -Loading MC -Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 into main simulation -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 -Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 now main simulation size is 420660 -Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 into main simulation -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 -Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 now main simulation size is 616803 -Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 into main simulation -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 -Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 now main simulation size is 740041 -Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 into main simulation -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 -Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 now main simulation size is 740161 -Loaded 740161 events in main simulation set -Loading Flux weighter -Loading nusquids objects -Weighting MC -Initializing simulation weights -Applying selfveto correction. -HESE reshuffle has not been performed. -Doing HESE reshuffle -Making data hist -Making sim hist -Constructing likelihood problem with default settings -Setting up Asimov set with parameters -Conventional normalization 1 -Prompt normalization 1 -Atmospheric muon normalization 1 -Astro component overall normalization 8 -Astro E component 0.333333 -Astro Mu component 0.333333 -Astro Tau component 0.333333 -Astroparticle balance 1 -Astro gamma 2.5 -Astro cutoff 10 -CR gamma 0 -Conv pi/k ratio 1 -Conv particle balance 1 -Conv zenith correction 0 -Dark normalization 0 -DOM efficiency 0.99 -Astro component overall normalization second 0 -Astro gamma second 2 -Holeice forward 0 -Anisotropy scale1 - -Remaking data hist -Reconstrucing likelihood problem -NULL min_llh 835.848062048 -NULL expectation [11.96478761 11.90169848 10.78840154 9.20158936 7.23605538 5.59952977 - 4.39850777 3.35653018 2.58491686 1.92834114 1.46813946 1.11581282 - 0.8492433 0.65154743 0.46727371 0.36211812 0.32718588 0.69746002 - 0.29475116 0.05150809] - -Applying new flavor physics to MC weights. -Entering ApplyNewPhysics -Entering ConstructNuSQuIDSObjectsForLV. -Rotated Unitary Transform -In for loop -Entering inner for loop -Leaving inner for loop diff --git a/test/NSI.out b/test/NSI.out deleted file mode 100644 index 657a9f3..0000000 --- a/test/NSI.out +++ /dev/null @@ -1,75 +0,0 @@ -test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 2u, std::allocator<double> > already registered; second conversion method ignored. - import GolemFitPy as gf -test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored. - import GolemFitPy as gf -reset_steering: 1 -reset_data: 1 -reset_npp: 1 -GolemFit constructor: checking paths -Constructing DiffuseWeightMaker -Loading data -Loading HESE data. -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/data/HESEData.txt -Loaded HESE 80 events. -Loading MC -Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 into main simulation -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 -Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 now main simulation size is 420660 -Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 into main simulation -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 -Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 now main simulation size is 616803 -Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 into main simulation -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 -Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 now main simulation size is 740041 -Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 into main simulation -Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 -Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 now main simulation size is 740161 -Loaded 740161 events in main simulation set -Loading Flux weighter -Loading nusquids objects -Weighting MC -Initializing simulation weights -Applying selfveto correction. -HESE reshuffle has not been performed. -Doing HESE reshuffle -Making data hist -Making sim hist -Constructing likelihood problem with default settings -Setting up Asimov set with parameters -Conventional normalization 1 -Prompt normalization 1 -Atmospheric muon normalization 1 -Astro component overall normalization 8 -Astro E component 0.333333 -Astro Mu component 0.333333 -Astro Tau component 0.333333 -Astroparticle balance 1 -Astro gamma 2.5 -Astro cutoff 10 -CR gamma 0 -Conv pi/k ratio 1 -Conv particle balance 1 -Conv zenith correction 0 -Dark normalization 0 -DOM efficiency 0.99 -Astro component overall normalization second 0 -Astro gamma second 2 -Holeice forward 0 -Anisotropy scale1 - -Remaking data hist -Reconstrucing likelihood problem -NULL min_llh 835.848062048 -NULL expectation [11.96478761 11.90169848 10.78840154 9.20158936 7.23605538 5.59952977 - 4.39850777 3.35653018 2.58491686 1.92834114 1.46813946 1.11581282 - 0.8492433 0.65154743 0.46727371 0.36211812 0.32718588 0.69746002 - 0.29475116 0.05150809] - -Applying new flavor physics to MC weights. -Applying New Physics in normal mode to NonStandardInteraction -HESE reshuffle has been performed. -0.1 mutau min_llh 836.352921667 -0.1 mutau expectation [16.63858821 16.04604012 14.28034835 11.80228198 9.06453127 6.8815312 - 5.32704912 4.02384379 3.07832186 2.28036523 1.73329938 1.30341591 - 0.99239868 0.75237435 0.53312845 0.4093463 0.35723267 0.71433954 - 0.30173276 0.05276177] diff --git a/test/test.png b/test/test.png Binary files differdeleted file mode 100644 index 8aac9e3..0000000 --- a/test/test.png +++ /dev/null diff --git a/test/test_LV.py b/test/test_LV.py deleted file mode 100644 index 96a1863..0000000 --- a/test/test_LV.py +++ /dev/null @@ -1,95 +0,0 @@ -#!/usr/bin/env python - -from __future__ import absolute_import, division - -import numpy as np -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt -from matplotlib import rc - -import GolemFitPy as gf - -rc('text', usetex=True) -rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) - -dp = gf.DataPaths() -steer = gf.SteeringParams() - -fig = plt.figure(figsize=[12, 8]) -ax = fig.add_subplot(111) -ax.set_xscale('log') -ax.set_yscale('log') - -npp = gf.NewPhysicsParams() -npp.type = gf.NewPhysicsType.None -# npp.n_lv = 1 -# npp.lambda_1 = 1.e100 -# npp.lambda_2 = 1.e100 - -golem = gf.GolemFit(dp, steer, npp) - -binning = golem.GetEnergyBinsMC() -ax.set_xlim(binning[0], binning[-1]) -# ax.set_ylim(binning[0], binning[-1]) - -fit_params = gf.FitParameters(gf.sampleTag.HESE) -golem.SetupAsimov(fit_params) - -exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) -ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, - drawstyle='steps-pre', label='NULL', linestyle='--') - -print 'NULL min_llh', golem.MinLLH().likelihood -print 'NULL expectation', exp -print - -npp = gf.NewPhysicsParams() -npp.type = gf.NewPhysicsType.LorentzViolation -npp.n_lv = 1 -npp.lambda_1 = 1.e20 -npp.lambda_2 = 1.e20 - -golem.SetNewPhysicsParams(npp) - -exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) -ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, - drawstyle='steps-pre', label='1e-20 LV', linestyle='--') - -print '1e20 LV min_llh', golem.MinLLH().likelihood -print '1e20 LV expectation', exp - -npp = gf.NewPhysicsParams() -npp.type = gf.NewPhysicsType.LorentzViolation -npp.n_lv = 1 -npp.lambda_1 = 1.e10 -npp.lambda_2 = 1.e10 - -golem.SetNewPhysicsParams(npp) - -exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) -ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, - drawstyle='steps-pre', label='1e-10 LV', linestyle='--') - -print '1e10 LV min_llh', golem.MinLLH().likelihood -print '1e10 LV expectation', exp - -npp = gf.NewPhysicsParams() -npp.type = gf.NewPhysicsType.LorentzViolation -npp.n_lv = 1 -npp.lambda_1 = 1.e-20 -npp.lambda_2 = 1.e-20 - -golem.SetNewPhysicsParams(npp) - -ax.tick_params(axis='x', labelsize=12) -ax.tick_params(axis='y', labelsize=12) -ax.set_xlabel(r'Deposited energy / GeV') -ax.set_ylabel(r'Events') -for xmaj in ax.xaxis.get_majorticklocs(): - ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1) -for ymaj in ax.yaxis.get_majorticklocs(): - ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1) - -legend = ax.legend(prop=dict(size=12)) -fig.savefig('test.png', bbox_inches='tight', dpi=250) diff --git a/test/test_NSI.png b/test/test_NSI.png Binary files differdeleted file mode 100644 index 000ea61..0000000 --- a/test/test_NSI.png +++ /dev/null diff --git a/test/test_NSI.py b/test/test_NSI.py deleted file mode 100644 index 9b49c93..0000000 --- a/test/test_NSI.py +++ /dev/null @@ -1,103 +0,0 @@ -#!/usr/bin/env python - -from __future__ import absolute_import, division - -import numpy as np -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt -from matplotlib import rc - -import GolemFitPy as gf - -rc('text', usetex=True) -rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) - -dp = gf.DataPaths() -steer = gf.SteeringParams() - -fig = plt.figure(figsize=[12, 8]) -ax = fig.add_subplot(111) -ax.set_xscale('log') -ax.set_yscale('log') - -npp = gf.NewPhysicsParams() -npp.type = gf.NewPhysicsType.None -# npp.epsilon_mutau = 0 -# npp.epsilon_prime = 0 - -golem = gf.GolemFit(dp, steer, npp) - -binning = golem.GetEnergyBinsMC() -ax.set_xlim(binning[0], binning[-1]) -# ax.set_ylim(binning[0], binning[-1]) - -fit_params = gf.FitParameters(gf.sampleTag.HESE) -golem.SetupAsimov(fit_params) - -exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) -ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, - drawstyle='steps-pre', label='NULL', linestyle='--') - -print 'NULL min_llh', golem.MinLLH().likelihood -print 'NULL expectation', exp -print - -npp = gf.NewPhysicsParams() -npp.type = gf.NewPhysicsType.NonStandardInteraction -npp.epsilon_mutau = 0.1 -# npp.epsilon_prime = 0 - -golem.SetNewPhysicsParams(npp) - -exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) -ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, - drawstyle='steps-pre', label='0.1 mutau', linestyle='--') - -print '0.1 mutau min_llh', golem.MinLLH().likelihood -print '0.1 mutau expectation', exp - -np.epsilon_mutau = 0.2 - -golem.SetNewPhysicsParams(npp) - -exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) -ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, - drawstyle='steps-pre', label='0.2 mutau', linestyle='--') - -print '0.2 mutau min_llh', golem.MinLLH().likelihood -print '0.2 mutau expectation', exp - -np.epsilon_mutau = 0.3 - -golem.SetNewPhysicsParams(npp) - -exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) -ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, - drawstyle='steps-pre', label='0.3 mutau', linestyle='--') - -print '0.3 mutau min_llh', golem.MinLLH().likelihood -print '0.3 mutau expectation', exp - -np.epsilon_mutau = 0.4 - -golem.SetNewPhysicsParams(npp) - -exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) -ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, - drawstyle='steps-pre', label='0.4 mutau', linestyle='--') - -print '0.4 mutau min_llh', golem.MinLLH().likelihood -print '0.4 mutau expectation', exp - -ax.tick_params(axis='x', labelsize=12) -ax.tick_params(axis='y', labelsize=12) -ax.set_xlabel(r'Deposited energy / GeV') -ax.set_ylabel(r'Events') -for xmaj in ax.xaxis.get_majorticklocs(): - ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1) -for ymaj in ax.yaxis.get_majorticklocs(): - ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1) - -legend = ax.legend(prop=dict(size=12)) -fig.savefig('test_NSI.png', bbox_inches='tight', dpi=250) diff --git a/utils/enums.py b/utils/enums.py index 90ca17d..359e104 100644 --- a/utils/enums.py +++ b/utils/enums.py @@ -47,6 +47,17 @@ class MCMCSeedType(Enum): GAUSSIAN = 2 +class SensitivityCateg(Enum): + FULL = 1 + FIXED_ANGLE = 2 + CORR_ANGLE = 3 + + +class StatCateg(Enum): + BAYESIAN = 1 + FREQUENTIST = 2 + + class SteeringCateg(Enum): P2_0 = 1 P2_1 = 2 diff --git a/utils/fr.py b/utils/fr.py index 8f71f78..6a405a7 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -9,11 +9,15 @@ Useful functions for the BSM flavour ratio analysis from __future__ import absolute_import, division -import sys +import argparse +from functools import partial import numpy as np from scipy import linalg +from utils.enums import EnergyDependance +from utils.misc import enum_parse, parse_bool + MASS_EIGENVALUES = [7.40E-23, 2.515E-21] """SM mass eigenvalues""" @@ -194,6 +198,83 @@ def normalise_fr(fr): return np.array(fr) / float(np.sum(fr)) +def estimate_scale(args, mass_eigenvalues=MASS_EIGENVALUES): + """Estimate the scale at which new physics will enter.""" + if args.energy_dependance is EnergyDependance.MONO: + scale = np.power( + 10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \ + np.log10(args.energy**(args.dimension-3)) + ) + elif args.energy_dependance is EnergyDependance.SPECTRAL: + scale = np.power( + 10, np.round( + np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \ + - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3)) + ) + ) + return scale + + +def fr_argparse(parser): + parser.add_argument( + '--measured-ratio', type=int, nargs=3, default=[1, 1, 1], + help='Set the central value for the measured flavour ratio at IceCube' + ) + parser.add_argument( + '--source-ratio', type=int, nargs=3, default=[2, 1, 0], + help='Set the source flavour ratio for the case when you want to fix it' + ) + parser.add_argument( + '--no-bsm', type=parse_bool, default='False', + help='Turn off BSM terms' + ) + parser.add_argument( + '--dimension', type=int, default=3, + help='Set the new physics dimension to consider' + ) + parser.add_argument( + '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance), + choices=EnergyDependance, + help='Type of energy dependance to use in the BSM fit' + ) + parser.add_argument( + '--spectral-index', default=-2, type=int, + help='Spectral index for spectral energy dependance' + ) + parser.add_argument( + '--binning', default=[1e4, 1e7, 5], type=float, nargs=3, + help='Binning for spectral energy dependance' + ) + parser.add_argument( + '--fix-source-ratio', type=parse_bool, default='False', + help='Fix the source flavour ratio' + ) + parser.add_argument( + '--fix-mixing', type=parse_bool, default='False', + help='Fix all mixing parameters to bi-maximal mixing' + ) + parser.add_argument( + '--fix-mixing-almost', type=parse_bool, default='False', + help='Fix all mixing parameters except s23' + ) + parser.add_argument( + '--fix-scale', type=parse_bool, default='False', + help='Fix the new physics scale' + ) + parser.add_argument( + '--scale', type=float, default=1e-30, + help='Set the new physics scale' + ) + parser.add_argument( + '--scale-region', type=float, default=1e10, + help='Set the size of the box to scan for the scale' + ) + parser.add_argument( + '--energy', type=float, default=1000, + help='Set the energy scale' + ) + + def fr_to_angles(ratios): """Convert from flavour ratio into the angular projection of the flavour ratios. @@ -218,7 +299,7 @@ NUFIT_U = angles_to_u((0.307, (1-0.02195)**2, 0.565, 3.97935)) def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES, nufit_u=NUFIT_U, no_bsm=False, fix_mixing=False, fix_mixing_almost=False, fix_scale=False, scale=None, - check_uni=True): + check_uni=True, epsilon=1e-9): """Construct the BSM mixing matrix from the BSM parameters. Parameters @@ -311,8 +392,8 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES, if check_uni: tu = test_unitarity(eg_vector) - if not abs(np.trace(tu) - 3.) < 1e-5 or \ - not abs(np.sum(tu) - 3.) < 1e-5: + if not abs(np.trace(tu) - 3.) < epsilon or \ + not abs(np.sum(tu) - 3.) < epsilon: raise AssertionError( 'Matrix is not unitary!\neg_vector\n{0}\ntest ' 'u\n{1}'.format(eg_vector, tu) @@ -378,4 +459,3 @@ def u_to_fr(source_fr, matrix): ) ratio = composition / np.sum(source_fr) return ratio - diff --git a/utils/gf.py b/utils/gf.py index 3f770eb..20afc75 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -9,7 +9,6 @@ Useful GolemFit wrappers for the BSM flavour ratio analysis from __future__ import absolute_import, division -import argparse import socket from functools import partial diff --git a/utils/likelihood.py b/utils/likelihood.py index c1208ab..1981c72 100644 --- a/utils/likelihood.py +++ b/utils/likelihood.py @@ -9,7 +9,6 @@ Likelihood functions for the BSM flavour ratio analysis from __future__ import absolute_import, division -import argparse from functools import partial import numpy as np @@ -34,6 +33,10 @@ def likelihood_argparse(parser): '--likelihood', default='gaussian', type=partial(enum_parse, c=Likelihood), choices=Likelihood, help='likelihood contour' ) + parser.add_argument( + '--sigma-ratio', type=float, default=0.01, + help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH' + ) def lnprior(theta, paramset): @@ -46,17 +49,17 @@ def lnprior(theta, paramset): return 0. -def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): +def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): """-Log likelihood function for a given theta.""" - if len(theta) != len(mcmc_paramset): + if len(theta) != len(llh_paramset): raise AssertionError( 'Length of MCMC scan is not the same as the input ' - 'params\ntheta={0}\nmcmc_paramset]{1}'.format(theta, mcmc_paramset) + 'params\ntheta={0}\nmcmc_paramset]{1}'.format(theta, llh_paramset) ) - for idx, param in enumerate(mcmc_paramset): + for idx, param in enumerate(llh_paramset): param.value = theta[idx] hypo_paramset = asimov_paramset - for param in mcmc_paramset.from_tag(ParamTag.NUISANCE): + for param in llh_paramset.from_tag(ParamTag.NUISANCE): hypo_paramset[param.name].value = param.value if args.energy_dependance is EnergyDependance.SPECTRAL: @@ -67,14 +70,14 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): if args.energy_dependance is EnergyDependance.MONO: source_flux = args.source_ratio elif args.energy_dependance is EnergyDependance.SPECTRAL: - source_flux = np.array( - [fr * np.power(bin_centers, args.spectral_index) - for fr in args.source_ratio] - ).T + source_flux = np.array( + [fr * np.power(bin_centers, args.spectral_index) + for fr in args.source_ratio] + ).T else: if args.energy_dependance is EnergyDependance.MONO: source_flux = fr_utils.angles_to_fr( - mcmc_paramset.from_tag(ParamTag.SRCANGLES, values=True) + llh_paramset.from_tag(ParamTag.SRCANGLES, values=True) ) elif args.energy_dependance is EnergyDependance.SPECTRAL: source_flux = np.array( @@ -82,7 +85,7 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): for fr in fr_utils.angles_to_fr(theta[-2:])] ).T - bsm_angles = mcmc_paramset.from_tag( + bsm_angles = llh_paramset.from_tag( [ParamTag.SCALE, ParamTag.MMANGLES], values=True ) @@ -134,11 +137,11 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): return gf_utils.get_llh_freq(fitter, hypo_paramset) -def ln_prob(theta, args, fitter, asimov_paramset, mcmc_paramset): - lp = lnprior(theta, paramset=mcmc_paramset) +def ln_prob(theta, args, fitter, asimov_paramset, llh_paramset): + lp = lnprior(theta, paramset=llh_paramset) if not np.isfinite(lp): return -np.inf return lp + triangle_llh( theta, args=args, asimov_paramset=asimov_paramset, - mcmc_paramset=mcmc_paramset, fitter=fitter + llh_paramset=llh_paramset, fitter=fitter ) diff --git a/utils/mcmc.py b/utils/mcmc.py index c712c3a..97b77dd 100644 --- a/utils/mcmc.py +++ b/utils/mcmc.py @@ -9,7 +9,6 @@ Useful functions to use an MCMC for the BSM flavour ratio analysis from __future__ import absolute_import, division -import argparse from functools import partial import emcee @@ -72,6 +71,14 @@ def mcmc_argparse(parser): type=partial(enum_parse, c=MCMCSeedType), choices=MCMCSeedType, help='Type of distrbution to make the initial MCMC seed' ) + parser.add_argument( + '--plot-angles', type=misc_utils.parse_bool, default='True', + help='Plot MCMC triangle in the angles space' + ) + parser.add_argument( + '--plot-elements', type=misc_utils.parse_bool, default='False', + help='Plot MCMC triangle in the mixing elements space' + ) def flat_seed(paramset, nwalkers): diff --git a/utils/misc.py b/utils/misc.py index f0c1ad4..59c5edb 100644 --- a/utils/misc.py +++ b/utils/misc.py @@ -9,7 +9,7 @@ Misc functions for the BSM flavour ratio analysis from __future__ import absolute_import, division -import os, sys +import os import errno import multiprocessing @@ -19,174 +19,7 @@ from operator import attrgetter import numpy as np -from utils.enums import Likelihood, ParamTag - - -class Param(object): - """Parameter class to store parameters. - """ - def __init__(self, name, value, ranges, seed=None, std=None, tex=None, tag=None): - self._seed = None - self._ranges = None - self._tex = None - self._tag = None - - self.name = name - self.value = value - self.seed = seed - self.ranges = ranges - self.std = std - self.tex = tex - self.tag = tag - - @property - def ranges(self): - return tuple(self._ranges) - - @ranges.setter - def ranges(self, values): - self._ranges = [val for val in values] - - @property - def seed(self): - if self._seed is None: return self.ranges - return tuple(self._seed) - - @seed.setter - def seed(self, values): - if values is None: return - self._seed = [val for val in values] - - @property - def tex(self): - return r'{0}'.format(self._tex) - - @tex.setter - def tex(self, t): - self._tex = t if t is not None else r'{\rm %s}' % self.name - - @property - def tag(self): - return self._tag - - @tag.setter - def tag(self, t): - if t is None: self._tag = ParamTag.NONE - else: - assert t in ParamTag - self._tag = t - - -class ParamSet(Sequence): - """Container class for a set of parameters. - """ - def __init__(self, *args): - param_sequence = [] - for arg in args: - try: - param_sequence.extend(arg) - except TypeError: - param_sequence.append(arg) - - if len(param_sequence) != 0: - # Disallow duplicated params - all_names = [p.name for p in param_sequence] - unique_names = set(all_names) - if len(unique_names) != len(all_names): - duplicates = set([x for x in all_names if all_names.count(x) > 1]) - raise ValueError('Duplicate definitions found for param(s): ' + - ', '.join(str(e) for e in duplicates)) - - # Elements of list must be Param type - assert all([isinstance(x, Param) for x in param_sequence]), \ - 'All params must be of type "Param"' - - self._params = param_sequence - - def __len__(self): - return len(self._params) - - def __getitem__(self, i): - if isinstance(i, int): - return self._params[i] - elif isinstance(i, basestring): - return self._by_name[i] - - def __getattr__(self, attr): - try: - return super(ParamSet, self).__getattribute__(attr) - except AttributeError: - t, v, tb = sys.exc_info() - try: - return self[attr] - except KeyError: - raise t, v, tb - - def __iter__(self): - return iter(self._params) - - def __str__(self): - o = '\n' - for obj in self._params: - o += '== {0:<15} = {1:<15}, tag={2:<15}\n'.format( - obj.name, obj.value, obj.tag - ) - return o - - @property - def _by_name(self): - return {obj.name: obj for obj in self._params} - - @property - def names(self): - return tuple([obj.name for obj in self._params]) - - @property - def labels(self): - return tuple([obj.tex for obj in self._params]) - - @property - def values(self): - return tuple([obj.value for obj in self._params]) - - @property - def seeds(self): - return tuple([obj.seed for obj in self._params]) - - @property - def ranges(self): - return tuple([obj.ranges for obj in self._params]) - - @property - def stds(self): - return tuple([obj.std for obj in self._params]) - - @property - def tags(self): - return tuple([obj.tag for obj in self._params]) - - @property - def params(self): - return self._params - - def to_dict(self): - return {obj.name: obj.value for obj in self._params} - - def from_tag(self, tag, values=False, index=False, invert=False): - if values and index: assert 0 - tag = np.atleast_1d(tag) - if not invert: - ps = [(idx, obj) for idx, obj in enumerate(self._params) - if obj.tag in tag] - else: - ps = [(idx, obj) for idx, obj in enumerate(self._params) - if obj.tag not in tag] - if values: - return tuple([io[1].value for io in ps]) - elif index: - return tuple([io[0] for io in ps]) - else: - return ParamSet([io[1] for io in ps]) +from utils.enums import Likelihood class SortingHelpFormatter(argparse.HelpFormatter): @@ -197,54 +30,32 @@ class SortingHelpFormatter(argparse.HelpFormatter): def gen_identifier(args): - mr = args.measured_ratio - si = args.sigma_ratio + f = '_DIM{0}'.format(args.dimension) + mr1, mr2, mr3 = args.measured_ratio if args.fix_source_ratio: - sr = args.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) + ) if args.fix_mixing: - out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fix_mixing'.format( - int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), - int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension - ) + f += '_fix_mixing' elif args.fix_mixing_almost: - out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fix_mixing_almost'.format( - int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), - int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension - ) + f += '_fix_mixing_almost' elif args.fix_scale: - out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fixed_scale_{8}'.format( - int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), - int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension, - args.scale - ) - else: - out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format( - int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000), - int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension - ) + f += '_fix_scale_{0}'.format(args.scale) else: - if args.fix_mixing: - out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fix_mixing'.format( - int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), - int(si*1000), args.dimension - ) - elif args.fix_mixing_almost: - out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fix_mixing_almost'.format( - int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), - int(si*1000), args.dimension - ) - elif args.fix_scale: - out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fixed_scale_{5}'.format( - int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), - int(si*1000), args.dimension, args.scale - ) - else: - out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}'.format( - int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), - int(si*1000), args.dimension - ) - if args.likelihood is Likelihood.FLAT: out += '_flat' - return out + f += '_mfr_{3:03d}_{4:03d}_{5:03d}'.format(mr1, mr2, mr3) + if args.fix_mixing: + f += '_fix_mixing' + elif args.fix_mixing_almost: + f += '_fix_mixing_almost' + elif args.fix_scale: + f += '_fix_scale_{0}'.format(args.scale) + if args.likelihood is Likelihood.FLAT: f += '_flat' + elif args.likelihood is Likelihood.GAUSSIAN: + f += '_sigma_{0:03d}'.format(int(args.sigma_ratio*1000)) + return f def gen_outfile_name(args): @@ -314,6 +125,13 @@ def make_dir(outfile): raise +def seed_parse(s): + if s.lower() == 'none': + return None + else: + return int(s) + + def thread_type(t): if t.lower() == 'max': return multiprocessing.cpu_count() diff --git a/utils/multinest.py b/utils/multinest.py new file mode 100644 index 0000000..3a7ab2d --- /dev/null +++ b/utils/multinest.py @@ -0,0 +1,117 @@ +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : April 19, 2018 + +""" +Useful functions to use MultiNest for the BSM flavour ratio analysis +""" + +from __future__ import absolute_import, division + +from functools import partial + +import numpy as np + +from pymultinest import analyse, run + +from utils.likelihood import triangle_llh +from utils.misc import make_dir, parse_bool + + +def CubePrior(cube, ndim, n_params): + # default are uniform priors + return ; + + +def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset, + args, fitter): + if ndim != len(mn_paramset): + raise AssertionError( + 'Length of MultiNest scan paramset is not the same as the input ' + 'params\ncube={0}\nmn_paramset]{1}'.format(cube, mn_paramset) + ) + pranges = mn_paramset.seeds + for i in xrange(ndim): + mn_paramset[i].value = (pranges[i][1]-pranges[i][0])*cube[i] + pranges[i][0] + for pm in mn_paramset.names: + llh_paramset[pm].value = mn_paramset[pm].value + theta = llh_paramset.values + # print 'llh_paramset', llh_paramset + return triangle_llh( + theta=theta, + args=args, + asimov_paramset=asimov_paramset, + llh_paramset=llh_paramset, + fitter=fitter + ) + + +def mn_argparse(parser): + parser.add_argument( + '--mn-run', type=parse_bool, default='True', + help='Run MultiNest' + ) + parser.add_argument( + '--mn-bins', type=int, default=10, + help='Number of bins for the Bayes factor plot' + ) + parser.add_argument( + '--mn-eval-bin', type=str, default='all', + help='Which bin to evalaute for Bayes factor plot' + ) + parser.add_argument( + '--mn-live-points', type=int, default=3000, + help='Number of live points for MultiNest evaluations' + ) + parser.add_argument( + '--mn-tolerance', type=float, default=0.01, + help='Tolerance for MultiNest' + ) + parser.add_argument( + '--mn-output', type=str, default='./mnrun/', + help='Folder to store MultiNest evaluations' + ) + parser.add_argument( + '--plot-multinest', type=parse_bool, default='False', + help='Plot MultiNest evidence' + ) + + +def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter): + """Run the MultiNest algorithm to calculate the evidence.""" + n_params = len(mn_paramset) + + for n in mn_paramset.names: + llh_paramset[n].value = mn_paramset[n].value + + lnProbEval = partial( + lnProb, + mn_paramset = mn_paramset, + llh_paramset = llh_paramset, + asimov_paramset = asimov_paramset, + args = args, + fitter = fitter + ) + + prefix = './mnrun/DIM{0}/{1:>019}_'.format( + args.dimension, np.random.randint(0, 2**63) + ) + make_dir(prefix) + print 'Running evidence calculation for {0}'.format(prefix) + result = run( + LogLikelihood = lnProbEval, + Prior = CubePrior, + n_dims = n_params, + n_live_points = args.mn_live_points, + evidence_tolerance = args.mn_tolerance, + outputfiles_basename = prefix, + importance_nested_sampling = True, + resume = False, + verbose = True + ) + + analyser = analyse.Analyzer( + outputfiles_basename=prefix, n_params=n_params + ) + return analyser.get_stats()['global evidence'] diff --git a/utils/param.py b/utils/param.py new file mode 100644 index 0000000..bf230df --- /dev/null +++ b/utils/param.py @@ -0,0 +1,235 @@ +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : April 19, 2018 + +""" +Param class and functions for the BSM flavour ratio analysis +""" + +from __future__ import absolute_import, division + +import sys + +from collections import Sequence + +import numpy as np + +from utils.plot import get_units +from utils.fr import fr_to_angles +from utils.enums import Likelihood, ParamTag + + +class Param(object): + """Parameter class to store parameters.""" + def __init__(self, name, value, ranges, seed=None, std=None, tex=None, tag=None): + self._seed = None + self._ranges = None + self._tex = None + self._tag = None + + self.name = name + self.value = value + self.seed = seed + self.ranges = ranges + self.std = std + self.tex = tex + self.tag = tag + + @property + def ranges(self): + return tuple(self._ranges) + + @ranges.setter + def ranges(self, values): + self._ranges = [val for val in values] + + @property + def seed(self): + if self._seed is None: return self.ranges + return tuple(self._seed) + + @seed.setter + def seed(self, values): + if values is None: return + self._seed = [val for val in values] + + @property + def tex(self): + return r'{0}'.format(self._tex) + + @tex.setter + def tex(self, t): + self._tex = t if t is not None else r'{\rm %s}' % self.name + + @property + def tag(self): + return self._tag + + @tag.setter + def tag(self, t): + if t is None: self._tag = ParamTag.NONE + else: + assert t in ParamTag + self._tag = t + + +class ParamSet(Sequence): + """Container class for a set of parameters.""" + def __init__(self, *args): + param_sequence = [] + for arg in args: + try: + param_sequence.extend(arg) + except TypeError: + param_sequence.append(arg) + + if len(param_sequence) != 0: + # Disallow duplicated params + all_names = [p.name for p in param_sequence] + unique_names = set(all_names) + if len(unique_names) != len(all_names): + duplicates = set([x for x in all_names if all_names.count(x) > 1]) + raise ValueError('Duplicate definitions found for param(s): ' + + ', '.join(str(e) for e in duplicates)) + + # Elements of list must be Param type + assert all([isinstance(x, Param) for x in param_sequence]), \ + 'All params must be of type "Param"' + + self._params = param_sequence + + def __len__(self): + return len(self._params) + + def __getitem__(self, i): + if isinstance(i, int): + return self._params[i] + elif isinstance(i, basestring): + return self._by_name[i] + + def __getattr__(self, attr): + try: + return super(ParamSet, self).__getattribute__(attr) + except AttributeError: + t, v, tb = sys.exc_info() + try: + return self[attr] + except KeyError: + raise t, v, tb + + def __iter__(self): + return iter(self._params) + + def __str__(self): + o = '\n' + for obj in self._params: + o += '== {0:<15} = {1:<15}, tag={2:<15}\n'.format( + obj.name, obj.value, obj.tag + ) + return o + + @property + def _by_name(self): + return {obj.name: obj for obj in self._params} + + @property + def names(self): + return tuple([obj.name for obj in self._params]) + + @property + def labels(self): + return tuple([obj.tex for obj in self._params]) + + @property + def values(self): + return tuple([obj.value for obj in self._params]) + + @property + def seeds(self): + return tuple([obj.seed for obj in self._params]) + + @property + def ranges(self): + return tuple([obj.ranges for obj in self._params]) + + @property + def stds(self): + return tuple([obj.std for obj in self._params]) + + @property + def tags(self): + return tuple([obj.tag for obj in self._params]) + + @property + def params(self): + return self._params + + def to_dict(self): + return {obj.name: obj.value for obj in self._params} + + def from_tag(self, tag, values=False, index=False, invert=False): + if values and index: assert 0 + tag = np.atleast_1d(tag) + if not invert: + ps = [(idx, obj) for idx, obj in enumerate(self._params) + if obj.tag in tag] + else: + ps = [(idx, obj) for idx, obj in enumerate(self._params) + if obj.tag not in tag] + if values: + return tuple([io[1].value for io in ps]) + elif index: + return tuple([io[0] for io in ps]) + else: + return ParamSet([io[1] for io in ps]) + + +def get_paramsets(args): + """Make the paramsets for generating the Asmimov MC sample and also running + the MCMC. + """ + asimov_paramset = [] + llh_paramset = [] + if args.likelihood == Likelihood.GOLEMFIT: + nuisance_paramset = define_nuisance() + asimov_paramset.extend(nuisance_paramset.params) + llh_paramset.extend(nuisance_paramset.params) + for parm in nuisance_paramset: + parm.value = args.__getattribute__(parm.name) + tag = ParamTag.BESTFIT + flavour_angles = fr_to_angles(args.measured_ratio) + asimov_paramset.extend([ + Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag), + Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag), + ]) + asimov_paramset = ParamSet(asimov_paramset) + + if not args.fix_mixing and not args.fix_mixing_almost: + tag = ParamTag.MMANGLES + e = 1e-10 + llh_paramset.extend([ + Param(name='s_12^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag), + Param(name='c_13^4', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag), + Param(name='s_23^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag), + Param(name='dcp', value=np.pi, ranges=[0.+e, 2*np.pi-e], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag) + ]) + if args.fix_mixing_almost: + tag = ParamTag.MMANGLES + llh_paramset.extend([ + Param(name='s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag) + ]) + 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 not args.fix_source_ratio: + tag = ParamTag.SRCANGLES + llh_paramset.extend([ + Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag), + Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag) + ]) + llh_paramset = ParamSet(llh_paramset) + return asimov_paramset, llh_paramset diff --git a/utils/plot.py b/utils/plot.py index 0b82675..66c888c 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -86,26 +86,6 @@ def interval(arr, percentile=68.): return sarr[curr_low], center, sarr[curr_up] -def plot_argparse(parser): - """Arguments for plotting.""" - parser.add_argument( - '--plot-angles', type=misc_utils.parse_bool, default='True', - help='Plot MCMC triangle in the angles space' - ) - parser.add_argument( - '--plot-elements', type=misc_utils.parse_bool, default='False', - help='Plot MCMC triangle in the mixing elements space' - ) - parser.add_argument( - '--plot-bayes', type=misc_utils.parse_bool, default='False', - help='Plot Bayes factor' - ) - parser.add_argument( - '--plot-angles-limit', type=misc_utils.parse_bool, default='False', - help='Plot limit vs BSM angles' - ) - - def flat_angles_to_u(x): """Convert from angles to mixing elements.""" return abs(angles_to_u(x)).astype(np.float32).flatten().tolist() @@ -138,32 +118,19 @@ def gen_figtext(args): mr1, mr2, mr3 = args.measured_ratio if args.fix_source_ratio: sr1, sr2, sr3 = args.source_ratio - if args.fix_scale: - t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \ - 'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \ - '{5:.2f}]\nDimension = {6}\nScale = {7}'.format( - sr1, sr2, sr3, mr1, mr2, mr3, args.dimension, - int(args.energy), args.scale - ) - else: - 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, mr1, mr2, mr3, args.dimension, - int(args.energy) - ) + 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, mr1, mr2, mr3, args.dimension, + int(args.energy) + ) else: - if args.fix_scale: - t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ - '{2:.2f}]\nDimension = {3}\nScale = {4}'.format( - mr1, mr2, mr3, args.dimension, int(args.energy), - args.scale - ) - else: - t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ - '{2:.2f}]\nDimension = {3}'.format( - mr1, mr2, mr3, args.dimension, int(args.energy) - ) + t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ + '{2:.2f}]\nDimension = {3}'.format( + mr1, mr2, mr3, args.dimension, int(args.energy) + ) + if args.fix_scale: + t += 'Scale = {0}'.format(args.scale) if args.likelihood is Likelihood.GAUSSIAN: t += '\nSigma = {0:.3f}'.format(args.sigma_ratio) if args.energy_dependance is EnergyDependance.SPECTRAL: @@ -176,7 +143,7 @@ def gen_figtext(args): return t -def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): +def chainer_plot(infile, outfile, outformat, args, llh_paramset): """Make the triangle plot.""" if not args.plot_angles and not args.plot_elements: return @@ -187,8 +154,8 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): misc_utils.make_dir(outfile) fig_text = gen_figtext(args) - axes_labels = mcmc_paramset.labels - ranges = mcmc_paramset.ranges + axes_labels = llh_paramset.labels + ranges = llh_paramset.ranges if args.plot_angles: print "Making triangle plots" @@ -208,7 +175,7 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): ), fontsize=10) # if not args.fix_mixing: - # sc_index = mcmc_paramset.from_tag(ParamTag.SCALE, index=True) + # 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 = ' @@ -222,11 +189,11 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): print "Making triangle plots" if args.fix_mixing_almost: raise NotImplementedError - nu_index = mcmc_paramset.from_tag(ParamTag.NUISANCE, index=True) - fr_index = mcmc_paramset.from_tag(ParamTag.MMANGLES, index=True) - sc_index = mcmc_paramset.from_tag(ParamTag.SCALE, index=True) + 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 = mcmc_paramset.from_tag(ParamTag.SRCANGLES, index=True) + 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])) @@ -264,45 +231,39 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): g.export(outfile+'_elements.'+of) -def bayes_factor_plot(dirname, outfile, outformat, args): - """Make Bayes factor plot.""" - if not args.plot_bayes: return - print "Making Bayes Factor plot" - print 'dirname', dirname +def plot_multinest(data, outfile, outformat, args, scale_param, label=None): + """Make MultiNest factor plot.""" + print "Making MultiNest Factor plot" fig_text = gen_figtext(args) + if label is not None: fig_text += '\n' + label - 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.sort(np.vstack(raw), axis=0) - print 'raw', raw - print 'raw.shape', raw.shape - scales, evidences = raw.T - null = evidences[0] + print 'data.shape', data.shape + scales, evidences = data.T + min_idx = np.argmin(scales) + null = evidences[min_idx] reduced_ev = -(evidences - null) fig = plt.figure(figsize=(7, 5)) ax = fig.add_subplot(111) - ax.set_xlim(np.log10(args.scale_region)) - ax.set_xlabel(r'${\rm log}_{10} \Lambda ' + get_units(args.dimension) +r'$') + ax.set_xlim(scale_param.ranges) + ax.set_xlabel('$'+scale_param.tex+'$') ax.set_ylabel(r'Bayes Factor') ax.plot(scales, reduced_ev) 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.3, 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.3, linewidth=1) at = AnchoredText( - fig_text, prop=dict(size=7), frameon=True, loc=2 + '\n'+fig_text, prop=dict(size=7), frameon=True, loc=2 ) - at.patch.set_boxstyle("round,pad=0.,rounding_size=0.5") + at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5") ax.add_artist(at) + misc_utils.make_dir(outfile) for of in outformat: fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) @@ -326,7 +287,8 @@ def plot_BSM_angles_limit(dirname, outfile, outformat, args, bayesian): for fn in filenames: if fn[-4:] == '.npy': raw.append(np.load(os.path.join(root, fn))) - raw = np.sort(np.vstack(raw), axis=0) + raw = np.vstack(raw) + raw = raw[np.argsort(raw[:,0])] print 'raw', raw print 'raw.shape', raw.shape sc_ranges = ( |
