diff options
Diffstat (limited to 'utils')
| -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 |
9 files changed, 541 insertions, 309 deletions
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 = ( |