diff options
Diffstat (limited to 'utils')
| -rw-r--r-- | utils/enums.py | 33 | ||||
| -rw-r--r-- | utils/fr.py | 140 | ||||
| -rw-r--r-- | utils/gf.py | 41 | ||||
| -rw-r--r-- | utils/llh.py | 104 | ||||
| -rw-r--r-- | utils/misc.py | 101 | ||||
| -rw-r--r-- | utils/mn.py | 13 | ||||
| -rw-r--r-- | utils/param.py | 66 | ||||
| -rw-r--r-- | utils/plot.py | 161 |
8 files changed, 203 insertions, 456 deletions
diff --git a/utils/enums.py b/utils/enums.py index 441a16f..e85158d 100644 --- a/utils/enums.py +++ b/utils/enums.py @@ -20,23 +20,9 @@ class DataType(Enum): REALISATION = 3 -class EnergyDependance(Enum): - MONO = 1 - SPECTRAL = 2 - - class Likelihood(Enum): - FLAT = 1 - GAUSSIAN = 2 - GOLEMFIT = 3 - GF_FREQ = 4 - - -class MixingScenario(Enum): - T12 = 1 - T13 = 2 - T23 = 3 - NONE = 4 + GOLEMFIT = 1 + GF_FREQ = 2 class ParamTag(Enum): @@ -60,14 +46,6 @@ class MCMCSeedType(Enum): GAUSSIAN = 2 -class SensitivityCateg(Enum): - FULL = 1 - FIXED_ANGLE = 2 - CORR_ANGLE = 3 - FIXED_ONE_ANGLE = 4 - CORR_ONE_ANGLE = 5 - - class StatCateg(Enum): BAYESIAN = 1 FREQUENTIST = 2 @@ -79,6 +57,7 @@ class SteeringCateg(Enum): class Texture(Enum): - OEU = 0 - OET = 1 - OUT = 2 + OEU = 1 + OET = 2 + OUT = 3 + NONE = 4 diff --git a/utils/fr.py b/utils/fr.py index 1d12fe5..b8eba44 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -13,7 +13,7 @@ from functools import partial import numpy as np -from utils.enums import EnergyDependance, MixingScenario +from utils.enums import Texture from utils.misc import enum_parse, parse_bool import mpmath as mp @@ -40,7 +40,17 @@ ASIN = mp.asin EXP = mp.exp MASS_EIGENVALUES = [7.40E-23, 2.515E-21] -"""SM mass eigenvalues""" +"""SM mass eigenvalues.""" + +SCALE_BOUNDARIES = { + 3 : (-32, -20), + 4 : (-40, -24), + 5 : (-48, -27), + 6 : (-56, -30), + 7 : (-64, -33), + 8 : (-72, -36) +} +"""Boundaries to scan the NP scale for each dimension.""" def determinant(x): @@ -244,38 +254,13 @@ def normalise_fr(fr): return np.array(fr) / float(np.sum(fr)) -def estimate_scale(args): - """Estimate the scale at which new physics will enter.""" - try: m_eign = args.m3x_2 - except: m_eign = MASS_EIGENVALUES[1] - if hasattr(args, 'scale'): - if args.scale != 0: - scale = args.scale - scale_region = (scale/args.scale_region, scale*args.scale_region) - return scale, scale_region - if args.energy_dependance is EnergyDependance.MONO: - scale = np.power( - 10, np.round(np.log10(m_eign/args.energy)) - \ - np.log10(args.energy**(args.dimension-3)) - ) - scale_region = (scale/args.scale_region, scale*args.scale_region) - elif args.energy_dependance is EnergyDependance.SPECTRAL: - lower_s = (m_eign/args.binning[-1]) / (args.binning[-1]**(args.dimension-3)) - upper_s = (m_eign/args.binning[0]) / (args.binning[0]**(args.dimension-3)) - scale = np.power(10, np.average(np.log10([lower_s, upper_s]))) - diff = upper_s / lower_s - scale_region = (lower_s/np.power(10, args.dimension), upper_s*diff*np.power(10, args.dimension)) - scale_region = [np.power(10, np.round(np.log10(x))) for x in scale_region] - return scale, scale_region - - def fr_argparse(parser): parser.add_argument( - '--measured-ratio', type=float, nargs=3, default=[1, 1, 1], - help='Set the central value for the measured flavour ratio at IceCube' + '--injected-ratio', type=float, nargs=3, required=False, + help='Injected ratio if not using data' ) parser.add_argument( - '--source-ratio', type=float, nargs=3, default=[2, 1, 0], + '--source-ratio', type=float, nargs=3, default=[1, 2, 0], help='Set the source flavour ratio for the case when you want to fix it' ) parser.add_argument( @@ -287,51 +272,13 @@ def fr_argparse(parser): 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=float, - help='Spectral index for spectral energy dependance' - ) - parser.add_argument( - '--fold-index', default='True', type=parse_bool, - help='Fold in the spectral index when using GolemFit' + '--texture', type=partial(enum_parse, c=Texture), + default='none', choices=Texture, help='Set the BSM mixing texture' ) parser.add_argument( '--binning', default=[6e4, 1e7, 20], 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=partial(enum_parse, c=MixingScenario), - default='None', choices=MixingScenario, - help='Fix all mixing parameters to choice of 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=0, - 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): @@ -363,8 +310,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, - sm_u=NUFIT_U, no_bsm=False, fix_mixing=MixingScenario.NONE, - fix_mixing_almost=False, fix_scale=False, scale=None, + sm_u=NUFIT_U, no_bsm=False, texture=Texture.NONE, check_uni=True, epsilon=1e-7): """Construct the BSM mixing matrix from the BSM parameters. @@ -388,17 +334,8 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES, no_bsm : bool Turn off BSM behaviour - fix_mixing : MixingScenario - Fix the BSM mixing angles - - fix_mixing_almost : bool - Fix the BSM mixing angles except one - - fix_scale : bool - Fix the BSM scale - - scale : float - Used with fix_scale - scale at which to fix + texture : Texture + BSM mixing texture check_uni : bool Check the resulting BSM mixing matrix is unitary @@ -422,33 +359,20 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES, 'got\n{0}'.format(sm_u) ) - if fix_mixing is not MixingScenario.NONE and fix_mixing_almost: - raise NotImplementedError( - '--fix-mixing and --fix-mixing-almost cannot be used together' - ) - if not isinstance(theta, (list, tuple)): theta = [theta] - if fix_mixing is MixingScenario.T12: - s12_2, c13_4, s23_2, dcp, sc2 = 0.5, 1.0, 0.0, 0., theta - elif fix_mixing is MixingScenario.T13: - s12_2, c13_4, s23_2, dcp, sc2 = 0.0, 0.25, 0.0, 0., theta - elif fix_mixing is MixingScenario.T23: - s12_2, c13_4, s23_2, dcp, sc2 = 0.0, 1.0, 0.5, 0., theta - elif fix_mixing_almost: - s12_2, c13_4, dcp = 0.5, 1.0-1E-6, 0. - s23_2, sc2 = theta - elif fix_scale: - s12_2, c13_4, s23_2, dcp = theta - sc2 = scale + z = 0.+1e-9 + if texture is Texture.OEU: + np_s12_2, np_c13_4, np_s23_2, np_dcp, sc2 = 0.5, 1.0, z, z, theta + elif texture is Texture.OET: + np_s12_2, np_c13_4, np_s23_2, np_dcp, sc2 = z, 0.25, z, z, theta + elif texture is Texture.OUT: + np_s12_2, np_c13_4, np_s23_2, np_dcp, sc2 = z, 1.0, 0.5, z, theta else: - s12_2, c13_4, s23_2, dcp, sc2 = theta + np_s12_2, np_c13_4, np_s23_2, np_dcp, sc2 = theta - if len(theta) != 0: - sc2 = np.power(10., sc2) - else: - sc2 = scale + sc2 = np.power(10., sc2) sc1 = sc2 / 100. mass_matrix = np.array( @@ -458,11 +382,11 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES, if no_bsm: eg_vector = cardano_eqn(sm_ham) else: - new_physics_u = angles_to_u((s12_2, c13_4, s23_2, dcp)) - scale_matrix = np.array( + NP_U = angles_to_u((np_s12_2, np_c13_4, np_s23_2, np_dcp)) + SC_U = np.array( [[0, 0, 0], [0, sc1, 0], [0, 0, sc2]] ) - bsm_term = (energy**(dim-3)) * np.dot(new_physics_u, np.dot(scale_matrix, new_physics_u.conj().T)) + bsm_term = (energy**(dim-3)) * np.dot(NP_U, np.dot(SC_U, NP_U.conj().T)) bsm_ham = sm_ham + bsm_term eg_vector = cardano_eqn(bsm_ham) diff --git a/utils/gf.py b/utils/gf.py index 13d5728..bc004bd 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -24,6 +24,9 @@ from utils.misc import enum_parse, parse_bool, thread_factors from utils.param import ParamSet +FITTER = None + + def fit_flags(llh_paramset): default_flags = { # False means it's not fixed in minimization @@ -74,9 +77,6 @@ def steering_params(args): elif args.likelihood is Likelihood.GF_FREQ: params.frequentist = True; - # For Tianlu - # params.years = [999] - if hasattr(args, 'binning'): params.minFitEnergy = args.binning[0] # GeV params.maxFitEnergy = args.binning[-1] # GeV @@ -95,59 +95,61 @@ def steering_params(args): return params -def setup_asimov(fitter, params): +def setup_asimov(params): print 'Injecting the model', params asimov_params = gf.FitParameters(gf.sampleTag.MagicTau) for parm in params: asimov_params.__setattr__(parm.name, float(parm.value)) - fitter.SetupAsimov(asimov_params) + FITTER.SetupAsimov(asimov_params) -def setup_realisation(fitter, params, seed): +def setup_realisation(params, seed): print 'Injecting the model', params asimov_params = gf.FitParameters(gf.sampleTag.MagicTau) for parm in params: asimov_params.__setattr__(parm.name, float(parm.value)) - fitter.Swallow(fitter.SpitRealization(asimov_params, seed)) + FITTER.Swallow(FITTER.SpitRealization(asimov_params, seed)) def setup_fitter(args, asimov_paramset): + global FITTER datapaths = gf.DataPaths() sparams = steering_params(args) npp = gf.NewPhysicsParams() - fitter = gf.GolemFit(datapaths, sparams, npp) + FITTER = gf.GolemFit(datapaths, sparams, npp) if args.data is DataType.ASIMOV: - setup_asimov(fitter, asimov_paramset) + setup_asimov(FITTER, asimov_paramset) elif args.data is DataType.REALISATION: seed = args.seed if args.seed is not None else 1 - setup_realisation(fitter, asimov_paramset, seed) + setup_realisation(FITTER, asimov_paramset, seed) elif args.data is DataType.REAL: print 'Using MagicTau DATA' - return fitter -def get_llh(fitter, params): +def get_llh(params): # print 'params', params fitparams = gf.FitParameters(gf.sampleTag.MagicTau) for parm in params: fitparams.__setattr__(parm.name, float(parm.value)) - llh = -fitter.EvalLLH(fitparams) + llh = -FITTER.EvalLLH(fitparams) return llh -def get_llh_freq(fitter, params): +def get_llh_freq(params): print 'setting to {0}'.format(params) fitparams = gf.FitParameters(gf.sampleTag.MagicTau) for parm in params: fitparams.__setattr__(parm.name, float(parm.value)) - fitter.SetFitParametersSeed(fitparams) - llh = -fitter.MinLLH().likelihood + FITTER.SetFitParametersSeed(fitparams) + llh = -FITTER.MinLLH().likelihood return llh -def data_distributions(fitter): - hdat = fitter.GetDataDistribution() - binedges = np.asarray([fitter.GetZenithBinsData(), fitter.GetEnergyBinsData()]) +def data_distributions(): + hdat = FITTER.GetDataDistribution() + binedges = np.asarray( + [FITTER.GetZenithBinsData(), FITTER.GetEnergyBinsData()] + ) return hdat, binedges @@ -164,4 +166,3 @@ def gf_argparse(parser): choices=SteeringCateg, help='use asimov/fake dataset with specific steering' ) - diff --git a/utils/llh.py b/utils/llh.py index 0c1e97d..93587b9 100644 --- a/utils/llh.py +++ b/utils/llh.py @@ -17,7 +17,7 @@ from scipy.stats import multivariate_normal, truncnorm from utils import fr as fr_utils from utils import gf as gf_utils -from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg +from utils.enums import Likelihood, ParamTag, PriorsCateg from utils.misc import enum_parse, gen_identifier, parse_bool @@ -34,22 +34,11 @@ def multi_gaussian(fr, fr_bf, sigma, offset=-320): return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr)) + offset -def likelihood_argparse(parser): +def llh_argparse(parser): parser.add_argument( - '--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' - ) - parser.add_argument( - '--save-measured-fr', type=parse_bool, default='False', - help='Output the measured flavour ratios' - ) - parser.add_argument( - '--output-measured-fr', type=str, default='./frs', - help='Output of the measured flavour ratios' + '--stat-method', default='bayesian', + type=partial(misc_utils.enum_parse, c=StatCateg), choices=StatCateg, + help='Statistical method to employ' ) @@ -82,7 +71,7 @@ def lnprior(theta, paramset): return prior -def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): +def triangle_llh(theta, args, asimov_paramset, llh_paramset): """Log likelihood function for a given theta.""" if len(theta) != len(llh_paramset): raise AssertionError( @@ -95,31 +84,14 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): for param in llh_paramset.from_tag(ParamTag.NUISANCE): hypo_paramset[param.name].value = param.value - if args.energy_dependance is EnergyDependance.SPECTRAL: - bin_centers = np.sqrt(args.binning[:-1]*args.binning[1:]) - bin_width = np.abs(np.diff(args.binning)) - if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ] \ - and args.fold_index: - args.spectral_index = -hypo_paramset['astroDeltaGamma'].value - - if args.fix_source_ratio: - 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 - else: - if args.energy_dependance is EnergyDependance.MONO: - source_flux = fr_utils.angles_to_fr( - llh_paramset.from_tag(ParamTag.SRCANGLES, values=True) - ) - elif args.energy_dependance is EnergyDependance.SPECTRAL: - source_flux = np.array( - [fr * np.power(bin_centers, args.spectral_index) - for fr in fr_utils.angles_to_fr(theta[-2:])] - ).T + bin_centers = np.sqrt(args.binning[:-1]*args.binning[1:]) + bin_width = np.abs(np.diff(args.binning)) + spectral_index = -hypo_paramset['astroDeltaGamma'].value + + source_flux = np.array( + [fr * np.power(bin_centers, spectral_index) + for fr in args.source_ratio] + ).T bsm_angles = llh_paramset.from_tag( [ParamTag.SCALE, ParamTag.MMANGLES], values=True @@ -139,21 +111,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): if args.no_bsm: fr = fr_utils.u_to_fr(source_flux, np.array(sm_u, dtype=np.complex256)) - elif args.energy_dependance is EnergyDependance.MONO: - u = fr_utils.params_to_BSMu( - theta = bsm_angles, - dim = args.dimension, - energy = args.energy, - mass_eigenvalues = mass_eigenvalues, - sm_u = sm_u, - no_bsm = args.no_bsm, - fix_mixing = args.fix_mixing, - fix_mixing_almost = args.fix_mixing_almost, - fix_scale = args.fix_scale, - scale = args.scale - ) - fr = fr_utils.u_to_fr(source_flux, u) - elif args.energy_dependance is EnergyDependance.SPECTRAL: + else: mf_perbin = [] for i_sf, sf_perbin in enumerate(source_flux): u = fr_utils.params_to_BSMu( @@ -163,10 +121,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): mass_eigenvalues = mass_eigenvalues, sm_u = sm_u, no_bsm = args.no_bsm, - fix_mixing = args.fix_mixing, - fix_mixing_almost = args.fix_mixing_almost, - fix_scale = args.fix_scale, - scale = args.scale + texture = args.texture, ) fr = fr_utils.u_to_fr(sf_perbin, u) mf_perbin.append(fr) @@ -181,35 +136,18 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): param.value = flavour_angles[idx] - if args.likelihood is Likelihood.FLAT: - llh = 1. - elif args.likelihood is Likelihood.GAUSSIAN: - fr_bf = args.measured_ratio - llh = multi_gaussian(fr, fr_bf, args.sigma_ratio) - elif args.likelihood is Likelihood.GOLEMFIT: - llh = gf_utils.get_llh(fitter, hypo_paramset) + if args.likelihood is Likelihood.GOLEMFIT: + llh = gf_utils.get_llh(hypo_paramset) elif args.likelihood is Likelihood.GF_FREQ: - llh = gf_utils.get_llh_freq(fitter, hypo_paramset) - - if args.save_measured_fr and args.burnin is False: - n = gen_identifier(args) + '.txt' - with open(args.output_measured_fr + n, 'a') as f: - f.write(r'{0:.3f} {1:.3f} {2:.3f}'.format( - float(fr[0]), float(fr[1]), float(fr[2]) - )) - for p in llh_paramset: - f.write(r' {0:.3f}'.format(p.value)) - f.write(' LLH = {0:.3f}'.format(llh)) - f.write('\n') - + llh = gf_utils.get_llh_freq(hypo_paramset) return llh -def ln_prob(theta, args, asimov_paramset, llh_paramset, fitter): +def ln_prob(theta, args, 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, - llh_paramset=llh_paramset, fitter=fitter + llh_paramset=llh_paramset ) diff --git a/utils/misc.py b/utils/misc.py index 26bd828..36d5330 100644 --- a/utils/misc.py +++ b/utils/misc.py @@ -19,7 +19,8 @@ from operator import attrgetter import numpy as np -from utils.enums import Likelihood, MixingScenario +from utils.enums import str_enum +from utils.enums import Likelihood, Texture class SortingHelpFormatter(argparse.HelpFormatter): @@ -31,30 +32,20 @@ class SortingHelpFormatter(argparse.HelpFormatter): def solve_ratio(fr): denominator = reduce(gcd, fr) - return [int(x/denominator) for x in fr] + f = [int(x/denominator) for x in fr] + if f[0] > 1E3 or f[1] > 1E3 or f[2] > 1E3: + return '{0:.2f}_{1:.2f}_{2:.2f}'.format(fr[0], fr[1], fr[2]) + else: + return '{0}_{1}_{2}'.format(f[0], f[1], f[2]) def gen_identifier(args): f = '_DIM{0}'.format(args.dimension) - mr1, mr2, mr3 = solve_ratio(args.measured_ratio) - if args.fix_source_ratio: - sr1, sr2, sr3 = solve_ratio(args.source_ratio) - f += '_sfr_{0:G}_{1:G}_{2:G}_mfr_{3:G}_{4:G}_{5:G}'.format( - sr1, sr2, sr3, mr1, mr2, mr3 - ) - else: - f += '_mfr_{3:03d}_{4:03d}_{5:03d}'.format(mr1, mr2, mr3) - - if args.fix_mixing is not MixingScenario.NONE: - f += '_{0}'.format(args.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)) + f += '_sfr_' + solve_ratio(args.source_ratio) + if args.data in [DataType.ASIMOV, DataType.REALISATION]: + f += '_mfr_' + solve_ratio(args.injected_ratio) + if args.Texture is not Texture.NONE: + f += '_{0}'.format(str_enum(args.texture)) return f @@ -163,3 +154,71 @@ def thread_factors(t): if t%x == 0: return (x, int(t/x)) + +def centers(x): + return (x[:-1]+x[1:])*0.5 + + +def get_units(dimension): + if dimension == 3: return r' / \:{\rm GeV}' + if dimension == 4: return r'' + if dimension == 5: return r' / \:{rm GeV}^{-1}' + if dimension == 6: return r' / \:{rm GeV}^{-2}' + if dimension == 7: return r' / \:{rm GeV}^{-3}' + if dimension == 8: return r' / \:{rm GeV}^{-4}' + + +def calc_nbins(x): + n = (np.max(x) - np.min(x)) / (2 * len(x)**(-1./3) * (np.percentile(x, 75) - np.percentile(x, 25))) + return np.floor(n) + + +def calc_bins(x): + nbins = calc_nbins(x) + return np.linspace(np.min(x), np.max(x)+2, num=nbins+1) + + +def most_likely(arr): + """Return the densest region given a 1D array of data.""" + binning = calc_bins(arr) + harr = np.histogram(arr, binning)[0] + return centers(binning)[np.argmax(harr)] + + +def interval(arr, percentile=68.): + """Returns the *percentile* shortest interval around the mode.""" + center = most_likely(arr) + sarr = sorted(arr) + delta = np.abs(sarr - center) + curr_low = np.argmin(delta) + curr_up = curr_low + npoints = len(sarr) + while curr_up - curr_low < percentile/100.*npoints: + if curr_low == 0: + curr_up += 1 + elif curr_up == npoints-1: + curr_low -= 1 + elif sarr[curr_up]-sarr[curr_low-1] < sarr[curr_up+1]-sarr[curr_low]: + curr_low -= 1 + elif sarr[curr_up]-sarr[curr_low-1] > sarr[curr_up+1]-sarr[curr_low]: + curr_up += 1 + elif (curr_up - curr_low) % 2: + # they are equal so step half of the time up and down + curr_low -= 1 + else: + curr_up += 1 + return sarr[curr_low], center, sarr[curr_up] + + +def flat_angles_to_u(x): + """Convert from angles to mixing elements.""" + return abs(angles_to_u(x)).astype(np.float32).flatten().tolist() + + +def myround(x, base=2, 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)) + + diff --git a/utils/mn.py b/utils/mn.py index 6582c80..c60d316 100644 --- a/utils/mn.py +++ b/utils/mn.py @@ -24,7 +24,7 @@ def CubePrior(cube, ndim, n_params): def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset, - args, fitter): + args): if ndim != len(mn_paramset): raise AssertionError( 'Length of MultiNest scan paramset is not the same as the input ' @@ -41,7 +41,6 @@ def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset, args=args, asimov_paramset=asimov_paramset, llh_paramset=llh_paramset, - fitter=fitter ) return llh @@ -61,13 +60,14 @@ def mn_argparse(parser): ) -def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter, +def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, identifier='mn'): """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 + print 'llh_paramset', llh_paramset lnProbEval = partial( lnProb, @@ -75,14 +75,13 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter, llh_paramset = llh_paramset, asimov_paramset = asimov_paramset, args = args, - fitter = fitter ) llh = '{0}'.format(args.likelihood).split('.')[1] data = '{0}'.format(args.data).split('.')[1] - sr1, sr2, sr3 = solve_ratio(args.source_ratio) - prefix = './mnrun/DIM{0}/{1}/{2}/s{3}{4}{5}/{6}'.format( - args.dimension, data, llh, sr1, sr2, sr3, identifier + src_string = solve_ratio(args.source_ratio) + prefix = './mnrun/DIM{0}/{1}/{2}/s{3}/{4}'.format( + args.dimension, data, llh, src_string, identifier ) make_dir(prefix) print 'Running evidence calculation for {0}'.format(prefix) diff --git a/utils/param.py b/utils/param.py index 572b65a..558018e 100644 --- a/utils/param.py +++ b/utils/param.py @@ -16,9 +16,8 @@ from copy import deepcopy import numpy as np -from utils.plot import get_units from utils.fr import fr_to_angles -from utils.enums import DataType, Likelihood, MixingScenario, ParamTag, PriorsCateg +from utils.enums import DataType, Likelihood, ParamTag, PriorsCateg class Param(object): @@ -219,66 +218,3 @@ class ParamSet(Sequence): elif isinstance(p, ParamSet): param_sequence.extend(p.params) return ParamSet(param_sequence) - - -def get_paramsets(args, nuisance_paramset): - """Make the paramsets for generating the Asmimov MC sample and also running - the MCMC. - """ - asimov_paramset = [] - llh_paramset = [] - - llh_paramset.extend( - [x for x in nuisance_paramset.from_tag(ParamTag.SM_ANGLES)] - ) - if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]: - gf_nuisance = [x for x in nuisance_paramset.from_tag(ParamTag.NUISANCE)] - asimov_paramset.extend(gf_nuisance) - llh_paramset.extend(gf_nuisance) - for parm in llh_paramset: - parm.value = args.__getattribute__(parm.name) - tag = ParamTag.BESTFIT - try: - flavour_angles = fr_to_angles(args.measured_ratio) - except: - flavour_angles = fr_to_angles(args.injected_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 hasattr(args, 'fix_source_ratio'): - if args.fix_mixing is MixingScenario.NONE and not args.fix_mixing_almost: - tag = ParamTag.MMANGLES - llh_paramset.extend([ - Param(name='np_s_12^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag), - Param(name='np_c_13^4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag), - Param(name='np_s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^2', tag=tag), - Param(name='np_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 - llh_paramset.extend([ - Param(name='np_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 - if hasattr(args, 'dimension'): - llh_paramset.append( - Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\left (\Lambda^{-1}'+get_units(args.dimension)+r'\right )', tag=tag) - ) - elif hasattr(args, 'dimensions'): - llh_paramset.append( - Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\left (\Lambda^{-1} / GeV^{-d+4}\right )', 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=[-1., 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 6161cfb..29ef5dc 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -19,12 +19,13 @@ from scipy.interpolate import splev, splprep from scipy.ndimage.filters import gaussian_filter import matplotlib as mpl +import matplotlib.patches as patches +import matplotlib.gridspec as gridspec mpl.use('Agg') + from matplotlib import rc from matplotlib import pyplot as plt from matplotlib.offsetbox import AnchoredText -import matplotlib.patches as patches -import matplotlib.gridspec as gridspec from matplotlib.lines import Line2D from matplotlib.patches import Patch @@ -32,22 +33,19 @@ from getdist import plots, mcsamples import ternary from ternary.heatmapping import polygon_generator -# print ternary.__file__ -# assert 0 import shapely.geometry as geometry from utils import misc as misc_utils from utils.enums import DataType, EnergyDependance, str_enum -from utils.enums import Likelihood, MixingScenario, ParamTag, StatCateg -from utils.enums import Texture +from utils.enums import Likelihood, ParamTag, StatCateg, Texture from utils.fr import angles_to_u, angles_to_fr -bayes_K = 1. # Substantial degree of belief. -# bayes_K = 3/2. # Strong degree of belief. -# bayes_K = 2. # Very strong degree of belief -# bayes_K = 5/2. # Decisive degree of belief +BAYES_K = 1. # Substantial degree of belief. +# BAYES_K = 3/2. # Strong degree of belief. +# BAYES_K = 2. # Very strong degree of belief +# BAYES_K = 5/2. # Decisive degree of belief if os.path.isfile('./plot_llh/paper.mplstyle'): @@ -58,64 +56,16 @@ if 'submitter' in socket.gethostname(): rc('text', usetex=False) -def centers(x): - return (x[:-1]+x[1:])*0.5 - - -def get_units(dimension): - if dimension == 3: return r' / \:{\rm GeV}' - if dimension == 4: return r'' - if dimension == 5: return r' / \:{rm GeV}^{-1}' - if dimension == 6: return r' / \:{rm GeV}^{-2}' - if dimension == 7: return r' / \:{rm GeV}^{-3}' - if dimension == 8: return r' / \:{rm GeV}^{-4}' - - -def calc_nbins(x): - n = (np.max(x) - np.min(x)) / (2 * len(x)**(-1./3) * (np.percentile(x, 75) - np.percentile(x, 25))) - return np.floor(n) - - -def calc_bins(x): - nbins = calc_nbins(x) - return np.linspace(np.min(x), np.max(x)+2, num=nbins+1) - - -def most_likely(arr): - """Return the densest region given a 1D array of data.""" - binning = calc_bins(arr) - harr = np.histogram(arr, binning)[0] - return centers(binning)[np.argmax(harr)] - - -def interval(arr, percentile=68.): - """Returns the *percentile* shortest interval around the mode.""" - center = most_likely(arr) - sarr = sorted(arr) - delta = np.abs(sarr - center) - curr_low = np.argmin(delta) - curr_up = curr_low - npoints = len(sarr) - while curr_up - curr_low < percentile/100.*npoints: - if curr_low == 0: - curr_up += 1 - elif curr_up == npoints-1: - curr_low -= 1 - elif sarr[curr_up]-sarr[curr_low-1] < sarr[curr_up+1]-sarr[curr_low]: - curr_low -= 1 - elif sarr[curr_up]-sarr[curr_low-1] > sarr[curr_up+1]-sarr[curr_low]: - curr_up += 1 - elif (curr_up - curr_low) % 2: - # they are equal so step half of the time up and down - curr_low -= 1 - else: - curr_up += 1 - return sarr[curr_low], center, sarr[curr_up] - - -def flat_angles_to_u(x): - """Convert from angles to mixing elements.""" - return abs(angles_to_u(x)).astype(np.float32).flatten().tolist() +def gen_figtext(args): + """Generate the figure text.""" + t = '' + t += 'Source flavour ratio = [{0}]'.format(solve_ratio(args.source_ratio)) + if args.data in [DataType.ASIMOV, DataType.REALISATION]: + t += '\nIC injected flavour ratio = [{0}]'.format( + solve_ratio(args.injected_ratio) + ) + t += '\nDimension = {0}'.format(args.dimension) + return t def plot_Tchain(Tchain, axes_labels, ranges): @@ -139,39 +89,6 @@ def plot_Tchain(Tchain, axes_labels, ranges): return g -def gen_figtext(args): - """Generate the figure text.""" - t = '' - mr1, mr2, mr3 = misc_utils.solve_ratio(args.measured_ratio) - if args.fix_source_ratio: - sr1, sr2, sr3 = misc_utils.solve_ratio(args.source_ratio) - t += 'Source flavour ratio = [{0}, {1}, {2}]'.format(sr1, sr2, sr3) - if args.data in [DataType.ASIMOV, DataType.REALISATION]: - t += '\nIC observed flavour ratio = [{0}, {1}, {2}]'.format( - mr1, mr2, mr3 - ) - t += '\nDimension = {0}'.format(args.dimension) - else: - if args.data in [DataType.ASIMOV, DataType.REALISATION]: - t += 'IC observed flavour ratio = [{0}, {1}, ' \ - '{2}]\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: - if not args.fold_index: - t += '\nSpectral Index = {0}'.format(int(args.spectral_index)) - t += '\nBinning = [{0}, {1}] TeV - {2} bins'.format( - int(args.binning[0]/1e3), int(args.binning[-1]/1e3), len(args.binning)-1 - ) - elif args.energy_dependance is EnergyDependance.MONO: - t += '\nEnergy = {0} TeV'.format(int(args.energy/1e3)) - return t - - def chainer_plot(infile, outfile, outformat, args, llh_paramset, fig_text=None): """Make the triangle plot.""" if hasattr(args, 'plot_elements'): @@ -287,13 +204,6 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset, fig_text=None): g.export(outfile+'_elements.'+of) -def myround(x, base=2, up=False, down=False): - if up == down and up is True: assert 0 - if up: return int(base * np.round(float(x)/base-0.5)) - elif down: return int(base * np.round(float(x)/base+0.5)) - else: int(base * np.round(float(x)/base)) - - def plot_statistic(data, outfile, outformat, args, scale_param, label=None): """Make MultiNest factor or LLH value plot.""" print 'Making Statistic plot' @@ -337,7 +247,7 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None): ax.plot(scales_rm, reduced_ev) - ax.axhline(y=np.log(10**(bayes_K)), color='red', alpha=1., linewidth=1.3) + ax.axhline(y=np.log(10**(BAYES_K)), color='red', alpha=1., linewidth=1.3) for ymaj in ax.yaxis.get_majorticklocs(): ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.3, linewidth=1) @@ -394,7 +304,7 @@ def plot_sens_full(data, outfile, outformat, args): null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: reduced_ev = -(statistic - null) - al = scales[reduced_ev > np.log(10**(bayes_K))] # Strong degree of belief + al = scales[reduced_ev > np.log(10**(BAYES_K))] # Strong degree of belief # al = scales[reduced_ev > 0.4] # Testing elif args.stat_method is StatCateg.FREQUENTIST: reduced_ev = -2*(statistic - null) @@ -437,8 +347,8 @@ def plot_sens_full(data, outfile, outformat, args): fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) -def plot_sens_fixed_angle_pretty(data, outfile, outformat, args): - print 'Making FIXED_ANGLE_PRETTY sensitivity plot' +def plot_table_sens(data, outfile, outformat, args): + print 'Making TABLE sensitivity plot' argsc = deepcopy(args) dims = len(data) LV_atmo_90pc_limits = { @@ -598,14 +508,14 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args): null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: reduced_ev = -(statistic_rm - null) - al = scales_rm[reduced_ev > np.log(10**(bayes_K))] + al = scales_rm[reduced_ev > np.log(10**(BAYES_K))] elif args.stat_method is StatCateg.FREQUENTIST: reduced_ev = -2*(statistic_rm - null) al = scales_rm[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks if len(al) == 0: print 'No points for DIM {0} FRS {1}!'.format(dim, src) continue - if reduced_ev[-1] < np.log(10**(bayes_K)) - 0.1: + if reduced_ev[-1] < np.log(10**(BAYES_K)) - 0.1: print 'Peaked contour does not exclude large scales! For ' \ 'DIM {0} FRS{1}!'.format(dim, src) continue @@ -739,14 +649,14 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: reduced_ev = -(statistic - null) - al = scales[reduced_ev > np.log(10**(bayes_K))] # Strong degree of belief + al = scales[reduced_ev > np.log(10**(BAYES_K))] # Strong degree of belief elif args.stat_method is StatCateg.FREQUENTIST: reduced_ev = -2*(statistic - null) al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks if len(al) == 0: print 'No points for DIM {0} FRS {1}!'.format(dim, src) continue - if reduced_ev[-1] < np.log(10**(bayes_K)) - 0.1: + if reduced_ev[-1] < np.log(10**(BAYES_K)) - 0.1: print 'Peaked contour does not exclude large scales! For ' \ 'DIM {0} FRS{1}!'.format(dim, src) continue @@ -903,7 +813,7 @@ def plot_sens_corr_angle(data, outfile, outformat, args): print sep_arrays if args.stat_method is StatCateg.BAYESIAN: - reduced_pdat_mask = (sep_arrays[2] > np.log(10**(bayes_K))) # Strong degree of belief + reduced_pdat_mask = (sep_arrays[2] > np.log(10**(BAYES_K))) # Strong degree of belief elif args.stat_method is StatCateg.FREQUENTIST: reduced_pdat_mask = (sep_arrays[2] > 4.61) # 90% CL for 2 DOFS via Wilks reduced_pdat = sep_arrays.T[reduced_pdat_mask].T @@ -1175,7 +1085,7 @@ def plot_source_ternary(data, outfile, outformat, args): for idim, dim in enumerate(args.dimensions): print '|||| DIM = {0}, {1}'.format(idim, dim) for isce in xrange(r_data.shape[1]): - print '|||| SCEN = {0}'.format(str_enum(Texture(isce))) + print '|||| SCEN = {0}'.format(str_enum(MixingScenario(isce+1))) fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(111) tax = get_tax(ax, scale=nsrcs) @@ -1204,14 +1114,14 @@ def plot_source_ternary(data, outfile, outformat, args): if args.stat_method is StatCateg.BAYESIAN: reduced_ev = -(statistic_rm - null) print 'reduced_ev', reduced_ev - al = scales_rm[reduced_ev > np.log(10**(bayes_K))] + al = scales_rm[reduced_ev > np.log(10**(BAYES_K))] else: assert 0 if len(al) == 0: print 'No points for DIM {0} FRS {1}!'.format(dim, src) interp_dict[src_sc] = -60 continue - if reduced_ev[-1] < np.log(10**(bayes_K)) - 0.1: + if reduced_ev[-1] < np.log(10**(BAYES_K)) - 0.1: print 'Peaked contour does not exclude large scales! For ' \ 'DIM {0} FRS{1}!'.format(dim, src) interp_dict[src_sc] = -60 @@ -1244,8 +1154,9 @@ def texture_label(x): raise AssertionError -def plot_source_ternary_1D(data, outfile, outformat, args): +def plot_x(data, outfile, outformat, args): """Ternary plot in the source flavour space for each operator texture.""" + print 'Making X sensitivity plot' sources = args.source_ratios x_1D = [] i_src_1D = [] @@ -1279,7 +1190,7 @@ def plot_source_ternary_1D(data, outfile, outformat, args): ax.set_ylabel(r'${\rm New\:Physics\:Scale}\:[\:{\rm log}_{10} (\Lambda_{d=6}^{-1}\:/\:{\rm GeV}^{-2})\: ]$', fontsize=18) ax.set_xlim(0, 1) for isce in xrange(r_data.shape[1]): - print '|||| SCEN = {0}'.format(str_enum(Texture(isce))) + print '|||| SCEN = {0}'.format(str_enum(MixingScenario(isce+1))) H = np.full(len(x_1D), np.nan) for ix, x in enumerate(x_1D): print '|||| X = {0}'.format(x) @@ -1303,13 +1214,13 @@ def plot_source_ternary_1D(data, outfile, outformat, args): if args.stat_method is StatCateg.BAYESIAN: reduced_ev = -(statistic_rm - null) print 'reduced_ev', reduced_ev - al = scales_rm[reduced_ev > np.log(10**(bayes_K))] + al = scales_rm[reduced_ev > np.log(10**(BAYES_K))] else: assert 0 if len(al) == 0: print 'No points for DIM {0} X {1}!'.format(dim, x) continue - if reduced_ev[-1] < np.log(10**(bayes_K)) - 0.1: + if reduced_ev[-1] < np.log(10**(BAYES_K)) - 0.1: print 'Peaked contour does not exclude large scales! For ' \ 'DIM {0} X {1}!'.format(dim, x) continue @@ -1320,7 +1231,7 @@ def plot_source_ternary_1D(data, outfile, outformat, args): H = ma.masked_invalid(H) H_0 = np.concatenate([[H[0]], H]) ax.step(be, H_0, linewidth=2, - label=texture_label(Texture(isce)), linestyle='-', + label=texture_label(MixingScenario(isce+1)), linestyle='-', drawstyle='steps-pre') print 'H', H print |