diff options
| -rwxr-xr-x | contour.py | 62 | ||||
| -rwxr-xr-x | fig2.py | 180 | ||||
| -rwxr-xr-x | fig2_austin.py | 162 | ||||
| -rwxr-xr-x | fr.py | 12 | ||||
| -rwxr-xr-x | plot_sens.py | 32 | ||||
| -rwxr-xr-x | sens.py | 28 | ||||
| -rw-r--r-- | submitter/sens_dag.py | 17 | ||||
| -rw-r--r-- | utils/enums.py | 6 | ||||
| -rw-r--r-- | utils/fr.py | 6 | ||||
| -rw-r--r-- | utils/gf.py | 5 | ||||
| -rw-r--r-- | utils/llh.py | 27 | ||||
| -rw-r--r-- | utils/plot.py | 197 |
12 files changed, 616 insertions, 118 deletions
@@ -10,6 +10,7 @@ HESE flavour ratio contour from __future__ import absolute_import, division +import os import argparse from functools import partial @@ -22,7 +23,7 @@ from utils import misc as misc_utils from utils import mn as mn_utils from utils import plot as plot_utils from utils.enums import str_enum -from utils.enums import DataType, Likelihood, ParamTag +from utils.enums import DataType, Likelihood, ParamTag, PriorsCateg from utils.param import Param, ParamSet, get_paramsets from pymultinest import Analyzer, run @@ -32,12 +33,13 @@ def define_nuisance(): """Define the nuisance parameters.""" nuisance = [] tag = ParamTag.NUISANCE + lg_prior = PriorsCateg.LIMITEDGAUSS nuisance.extend([ - 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=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0.1, 10.], std=0.4, prior=lg_prior, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 20.], std=2.4, prior=lg_prior, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 10.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0., 5. ], ranges=[0. , 20.], std=1.5, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) ]) return ParamSet(nuisance) @@ -120,12 +122,12 @@ def gen_identifier(args): def gen_figtext(args, asimov_paramset): f = '' if args.data is DataType.REAL: - f += 'IceCube Preliminary\n' + f += 'IceCube Preliminary' else: ir1, ir2, ir3 = misc_utils.solve_ratio(args.injected_ratio) - f += 'Injected ratio = [{0}, {1}, {2}]\n'.format(ir1, ir2, ir3) + f += 'Injected ratio = [{0}, {1}, {2}]'.format(ir1, ir2, ir3) for param in asimov_paramset: - f += 'Injected {0} = {1:.3f}\n'.format( + f += '\nInjected {0:20s} = {1:.3f}'.format( param.name, param.nominal_value ) return f @@ -210,6 +212,9 @@ def main(): fitter = fitter ) + cwd = os.getcwd() + os.chdir(prefix[:-len(os.path.basename(prefix))]) + print 'Running evidence calculation for {0}'.format(prefix) run( LogLikelihood = lnProbEval, @@ -217,12 +222,14 @@ def main(): n_dims = n_params, n_live_points = args.mn_live_points, evidence_tolerance = args.mn_tolerance, - outputfiles_basename = prefix, + outputfiles_basename = prefix[-len(os.path.basename(prefix)):], importance_nested_sampling = True, resume = False, verbose = True ) + os.chdir(cwd) + # Analyze analyser = Analyzer( outputfiles_basename=prefix, n_params=n_params @@ -230,38 +237,53 @@ def main(): print analyser pranges = hypo_paramset.ranges + + bf = analyser.get_best_fit()['parameters'] + for i in xrange(len(bf)): + bf[i] = (pranges[i][1]-pranges[i][0])*bf[i] + pranges[i][0] + print 'bestfit = ', bf + print 'bestfit log_likelihood', analyser.get_best_fit()['log_likelihood'] + + print + print '{0:50} = {1}'.format('global evidence', analyser.get_stats()['global evidence']) + print + fig_text = gen_figtext(args, asimov_paramset) + fig_text += '\nBestfit LLH = {0}'.format(analyser.get_best_fit()['log_likelihood']) + fig_text += '\nBestfits = ' + for x in bf: fig_text += '{0:.2f} '.format(x) - if args.plot_chains: + if args.plot_chains or args.plot_triangle: chains = analyser.get_data()[:,2:] for x in chains: for i in xrange(len(x)): x[i] = (pranges[i][1]-pranges[i][0])*x[i] + pranges[i][0] + + if args.plot_chains: + of = outfile[:5]+outfile[5:].replace('data', 'plots')+'_posterior' plot_utils.chainer_plot( infile = chains, - outfile = outfile[:5]+outfile[5:].replace('data', 'plots')+'_posterior', - outformat = ['pdf'], + outfile = of, + outformat = ['png'], args = args, llh_paramset = hypo_paramset, fig_text = fig_text ) + print 'Saved plot', of if args.plot_triangle: - llh = analyser.get_data()[:,1] + llh = -0.5 * analyser.get_data()[:,1] - chains = analyser.get_data()[:,2:] - for x in chains: - for i in xrange(len(x)): - x[i] = (pranges[i][1]-pranges[i][0])*x[i] + pranges[i][0] flavour_angles = chains[:,-2:] flavour_ratios = np.array( - map(fr_utils.angles_to_fr, flavour_angles), dtype=np.float + map(fr_utils.angles_to_fr, flavour_angles) ) + of = outfile[:5]+outfile[5:].replace('data', 'plots')+'_triangle' plot_utils.triangle_project( frs = flavour_ratios, llh = llh, - outfile = outfile[:5]+outfile[5:].replace('data', 'plots')+'_triangle', + outfile = of, outformat = ['png'], args = args, llh_paramset = hypo_paramset, @@ -0,0 +1,180 @@ +#! /usr/bin/env python +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : February 24, 2019 + +""" +HESE BSM Flavour Figure 2 +""" + +from __future__ import absolute_import, division + +import argparse +from functools import partial + +import numpy as np + +from utils import fr as fr_utils +from utils import misc as misc_utils +from utils import plot as plot_utils +from utils.enums import str_enum +from utils.enums import Likelihood, ParamTag, PriorsCateg +from utils.param import Param, ParamSet + +from matplotlib import pyplot as plt + +from pymultinest import Analyzer + + +def define_nuisance(): + """Define the nuisance parameters.""" + nuisance = [] + tag = ParamTag.NUISANCE + lg_prior = PriorsCateg.LIMITEDGAUSS + nuisance.extend([ + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0.1, 10.], std=0.4, prior=lg_prior, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 20.], std=2.4, prior=lg_prior, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 10.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0., 5. ], ranges=[0. , 20.], std=1.5, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + ]) + return ParamSet(nuisance) + + +def get_paramsets(args, nuisance_paramset): + paramset = [] + if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]: + gf_nuisance = [x for x in nuisance_paramset.from_tag(ParamTag.NUISANCE)] + paramset.extend(gf_nuisance) + tag = ParamTag.BESTFIT + paramset.extend([ + Param(name='astroFlavorAngle1', value=0, ranges=[0., 1.], std=0.2, tag=tag), + Param(name='astroFlavorAngle2', value=0, ranges=[-1., 1.], std=0.2, tag=tag), + ]) + paramset = ParamSet(paramset) + return paramset + + +def process_args(args): + """Process the input args.""" + if args.likelihood is not Likelihood.GOLEMFIT \ + and args.likelihood is not Likelihood.GF_FREQ: + raise AssertionError( + 'Likelihood method {0} not supported for this ' + 'script!\nChoose either GOLEMFIT or GF_FREQ'.format( + str_enum(args.likelihood) + ) + ) + +def parse_args(args=None): + """Parse command line arguments""" + parser = argparse.ArgumentParser( + description="HESE BSM Flavour Figure 2", + formatter_class=misc_utils.SortingHelpFormatter, + ) + parser.add_argument( + '--likelihood', default='golemfit', + type=partial(misc_utils.enum_parse, c=Likelihood), + choices=Likelihood, help='likelihood contour' + ) + parser.add_argument( + '--contour-dir', type=str, + help='Path to directory containing MultiNest runs' + ) + parser.add_argument( + '--outfile', type=str, default='./untitled', + help='Output path' + ) + 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) + + paramset = get_paramsets(args, define_nuisance()) + n_params = len(paramset) + print n_params + + # Data + data_path = '{0}/{1}/real'.format( + args.contour_dir, str_enum(args.likelihood).lower() + ) + prefix = '{0}/_{1}_REAL_mn_'.format(data_path, str_enum(args.likelihood)) + analyser = Analyzer( + outputfiles_basename=prefix, n_params=n_params + ) + print analyser + + pranges = paramset.ranges + + bf = analyser.get_best_fit()['parameters'] + for i in xrange(len(bf)): + bf[i] = (pranges[i][1]-pranges[i][0])*bf[i] + pranges[i][0] + print 'bestfit = ', bf + print 'bestfit log_likelihood', analyser.get_best_fit()['log_likelihood'] + + print + print '{0:50} = {1}'.format('global evidence', analyser.get_stats()['global evidence']) + print + + chains = analyser.get_data()[:,2:] + for x in chains: + for i in xrange(len(x)): + x[i] = (pranges[i][1]-pranges[i][0])*x[i] + pranges[i][0] + + llh = -0.5 * analyser.get_data()[:,1] + + flavour_angles = chains[:,-2:] + flavour_ratios = np.array( + map(fr_utils.angles_to_fr, flavour_angles) + ) + + nbins = 25 + + fig = plt.figure(figsize=(8, 8)) + ax = fig.add_subplot(111) + tax = plot_utils.get_tax(ax, scale=nbins) + + plot_utils.flavour_contour( + frs = flavour_ratios, + ax = ax, + nbins = nbins, + coverage = 99, + linewidth = 2, + color = 'green' + ) + + plot_utils.flavour_contour( + frs = flavour_ratios, + ax = ax, + nbins = nbins, + coverage = 90, + linewidth = 2, + color = 'blue' + ) + + plot_utils.flavour_contour( + frs = flavour_ratios, + ax = ax, + nbins = nbins, + coverage = 68, + linewidth = 2, + color = 'red' + ) + + ax.legend() + + fig.savefig('test.png', bbox_inches='tight', dpi=150) + + print "DONE!" + + +main.__doc__ = __doc__ + + +if __name__ == '__main__': + main() diff --git a/fig2_austin.py b/fig2_austin.py new file mode 100755 index 0000000..8a4cc01 --- /dev/null +++ b/fig2_austin.py @@ -0,0 +1,162 @@ +#! /usr/bin/env python +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : February 24, 2019 + +""" +HESE BSM Flavour Figure 2 +""" + +from __future__ import absolute_import, division + +import argparse +from functools import partial + +import numpy as np + +from utils import fr as fr_utils +from utils import misc as misc_utils +from utils import plot as plot_utils +from utils.enums import str_enum +from utils.enums import Likelihood, ParamTag, PriorsCateg +from utils.param import Param, ParamSet + +from matplotlib import pyplot as plt + +# from pymultinest import Analyzer +import json + + +def define_nuisance(): + """Define the nuisance parameters.""" + nuisance = [] + tag = ParamTag.NUISANCE + lg_prior = PriorsCateg.LIMITEDGAUSS + nuisance.extend([ + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0.1, 10.], std=0.4, prior=lg_prior, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 20.], std=2.4, prior=lg_prior, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 10.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0., 5. ], ranges=[0. , 20.], std=1.5, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + ]) + return ParamSet(nuisance) + + +def get_paramsets(args, nuisance_paramset): + paramset = [] + if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]: + gf_nuisance = [x for x in nuisance_paramset.from_tag(ParamTag.NUISANCE)] + paramset.extend(gf_nuisance) + tag = ParamTag.BESTFIT + paramset.extend([ + Param(name='astroFlavorAngle1', value=0, ranges=[0., 1.], std=0.2, tag=tag), + Param(name='astroFlavorAngle2', value=0, ranges=[-1., 1.], std=0.2, tag=tag), + ]) + paramset = ParamSet(paramset) + return paramset + + +def process_args(args): + """Process the input args.""" + if args.likelihood is not Likelihood.GOLEMFIT \ + and args.likelihood is not Likelihood.GF_FREQ: + raise AssertionError( + 'Likelihood method {0} not supported for this ' + 'script!\nChoose either GOLEMFIT or GF_FREQ'.format( + str_enum(args.likelihood) + ) + ) + +def parse_args(args=None): + """Parse command line arguments""" + parser = argparse.ArgumentParser( + description="HESE BSM Flavour Figure 2", + formatter_class=misc_utils.SortingHelpFormatter, + ) + parser.add_argument( + '--likelihood', default='golemfit', + type=partial(misc_utils.enum_parse, c=Likelihood), + choices=Likelihood, help='likelihood contour' + ) + parser.add_argument( + '--contour-dir', type=str, + help='Path to directory containing MultiNest runs' + ) + parser.add_argument( + '--outfile', type=str, default='./untitled', + help='Output path' + ) + 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) + + paramset = get_paramsets(args, define_nuisance()) + n_params = len(paramset) + print n_params + + # Data + data_path = '/home/aschneider/programs/GOLEMSPACE/sources/GolemFit/scripts/diffuse/mcmcs/results/dpl_numu_prior_flavor_20190302-162221-a747f528-8aa6-4488-8c80-059572c099fe.json' + with open(data_path) as f: + d_json = json.load(f) + + names = d_json['func_args'] + chains = np.array(d_json['chain']) + print 'names', names + print 'chains.shape', chains.shape + + flavour_angles = chains[:,5:7] + flavour_ratios = np.array( + map(fr_utils.angles_to_fr, flavour_angles) + ) + + nbins = 25 + + fig = plt.figure(figsize=(8, 8)) + ax = fig.add_subplot(111) + tax = plot_utils.get_tax(ax, scale=nbins) + + plot_utils.flavour_contour( + frs = flavour_ratios, + ax = ax, + nbins = nbins, + coverage = 99, + linewidth = 2, + color = 'green' + ) + + plot_utils.flavour_contour( + frs = flavour_ratios, + ax = ax, + nbins = nbins, + coverage = 90, + linewidth = 2, + color = 'blue' + ) + + plot_utils.flavour_contour( + frs = flavour_ratios, + ax = ax, + nbins = nbins, + coverage = 68, + linewidth = 2, + color = 'red' + ) + + ax.legend() + + fig.savefig('test_austin.png', bbox_inches='tight', dpi=150) + + print "DONE!" + + +main.__doc__ = __doc__ + + +if __name__ == '__main__': + main() @@ -32,7 +32,7 @@ def define_nuisance(): tag = ParamTag.SM_ANGLES nuisance = [] g_prior = PriorsCateg.GAUSSIAN - hg_prior = PriorsCateg.HALFGAUSS + hg_prior = PriorsCateg.LIMITEDGAUSS e = 1e-9 nuisance.extend([ Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), @@ -62,11 +62,11 @@ def define_nuisance(): ]) tag = ParamTag.NUISANCE nuisance.extend([ - 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) + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0.1, 10.], std=0.4, prior=lg_prior, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 20.], std=2.4, prior=lg_prior, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 10.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0., 5. ], ranges=[0. , 20.], std=1.5, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) ]) return ParamSet(nuisance) diff --git a/plot_sens.py b/plot_sens.py index 2e3bcd0..b12389f 100755 --- a/plot_sens.py +++ b/plot_sens.py @@ -34,13 +34,13 @@ def define_nuisance(): """Define the nuisance parameters.""" tag = ParamTag.SM_ANGLES g_prior = PriorsCateg.GAUSSIAN - hg_prior = PriorsCateg.HALFGAUSS + lg_prior = PriorsCateg.LIMITEDGAUSS e = 1e-9 nuisance = [ - Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), - Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.995, 1-e], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=hg_prior, tag=tag), - Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), - Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), + Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=lg_prior, tag=tag), + Param(name='c_13_4', value=(1-(0.02206))**2, seed=[0.950, 0.961], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=lg_prior, tag=tag), + Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=lg_prior, tag=tag), + Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), Param( name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23], std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag @@ -52,11 +52,11 @@ def define_nuisance(): ] tag = ParamTag.NUISANCE nuisance.extend([ - 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) + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0.1, 10.], std=0.4, prior=lg_prior, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 20.], std=2.4, prior=lg_prior, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 10.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0., 5. ], ranges=[0. , 20.], std=1.5, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) ]) return ParamSet(nuisance) @@ -233,11 +233,13 @@ def main(): infile += '/gaussian/' if args.likelihood is Likelihood.GAUSSIAN: infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_')) - # infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/fr_stat'.format( - infile += '/DIM{0}/fix_ifr/prior/{1}/{2}/{3}/fr_stat'.format( + infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/fr_stat'.format( + # infile += '/DIM{0}/fix_ifr/prior/{1}/{2}/{3}/fr_stat'.format( # infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/old/fr_stat'.format( # infile += '/DIM{0}/fix_ifr/seed2/{1}/{2}/{3}/fr_stat'.format( # infile += '/DIM{0}/fix_ifr/100TeV/{1}/{2}/{3}/fr_stat'.format( + # infile += '/DIM{0}/fix_ifr/strictprior/{1}/{2}/{3}/fr_stat'.format( + # infile += '/DIM{0}/fix_ifr/noprior/{1}/{2}/{3}/fr_stat'.format( dim, *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data]) ) + misc_utils.gen_identifier(argsc) print '== {0:<25} = {1}'.format('infile', infile) @@ -283,10 +285,12 @@ def main(): base_infile += '/gaussian/' if args.likelihood is Likelihood.GAUSSIAN: base_infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_')) - # base_infile += '/DIM{0}/fix_ifr'.format(dim) - base_infile += '/DIM{0}/fix_ifr/prior'.format(dim) + base_infile += '/DIM{0}/fix_ifr'.format(dim) + # base_infile += '/DIM{0}/fix_ifr/prior'.format(dim) # base_infile += '/DIM{0}/fix_ifr/seed2'.format(dim) # base_infile += '/DIM{0}/fix_ifr/100TeV'.format(dim) + # base_infile += '/DIM{0}/fix_ifr/strictprior'.format(dim) + # base_infile += '/DIM{0}/fix_ifr/noprior'.format(dim) for isrc, src in enumerate(args.source_ratios): argsc.source_ratio = src @@ -35,11 +35,12 @@ def define_nuisance(): tag = ParamTag.SM_ANGLES nuisance = [] g_prior = PriorsCateg.GAUSSIAN + lg_prior = PriorsCateg.LIMITEDGAUSS e = 1e-9 nuisance.extend([ - Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), - Param(name='c_13_4', value=(1-(0.02206))**2, seed=[0.950, 0.961], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=g_prior, tag=tag), - Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), + Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=lg_prior, tag=tag), + Param(name='c_13_4', value=(1-(0.02206))**2, seed=[0.950, 0.961], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=lg_prior, tag=tag), + Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=lg_prior, tag=tag), Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), Param( name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23], @@ -52,11 +53,21 @@ def define_nuisance(): ]) tag = ParamTag.NUISANCE nuisance.extend([ - 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=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0.1, 10.], std=0.4, prior=lg_prior, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 20.], std=2.4, prior=lg_prior, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 10.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0., 5. ], ranges=[0. , 20.], std=1.5, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + # 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=[2.4, 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + # Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0.5, 2. ], std=0.3, tag=tag), + # Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 6. ], std=0.05, tag=tag), + # Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0.1, 2. ], std=0.1, tag=tag), + # Param(name='astroNorm', value=6.9, seed=[0.1, 10.], ranges=[0.1, 10.], std=0.1, tag=tag), + # Param(name='astroDeltaGamma', value=2.5, seed=[2.4, 3. ], ranges=[2.4, 3. ], std=0.1, tag=tag) ]) return ParamSet(nuisance) @@ -299,6 +310,7 @@ def main(): continue print '## Evidence = {0}'.format(stat) elif args.stat_method is StatCateg.FREQUENTIST: + raise NotImplementedError('Still needs testing') def fn(x): # TODO(shivesh): should be seed or ranges? # Force prior ranges to be inside "seed" diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py index bfab194..7423d34 100644 --- a/submitter/sens_dag.py +++ b/submitter/sens_dag.py @@ -35,17 +35,18 @@ GLOBAL_PARAMS.update(dict( # MultiNest GLOBAL_PARAMS.update(dict( - # mn_live_points = 1000, - mn_live_points = 500, + mn_live_points = 1000, + # mn_live_points = 500, + # mn_live_points = 300, # mn_tolerance = 0.1, mn_tolerance = 0.3, mn_output = './mnrun' )) # FR -dimension = [6] +# dimension = [6] # dimension = [3, 6] -# dimension = [3, 4, 5, 6, 7, 8] +dimension = [3, 4, 5, 6, 7, 8] GLOBAL_PARAMS.update(dict( threads = 1, binning = '6e4 1e7 20', @@ -84,7 +85,9 @@ outfile = 'dagman_FR_SENS_{0}_{1}_{2}_{3}'.format( # outfile += '_seed2' # outfile += '_tol03' # outfile += '_NULL' -outfile += '_prior' +# outfile += '_prior' +# outfile += '_strictprior' +# outfile += '_noprior' outfile += '.submit' golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit' condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub' @@ -112,7 +115,9 @@ with open(outfile, 'w') as f: # output += 'seed2/' # output += 'mn_noverlap/' # output += 'tol_03/' - output += 'prior/' + # output += 'prior/' + # output += 'strictprior/' + # output += 'noprior/' for r in xrange(sens_runs): print 'run', r f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script)) diff --git a/utils/enums.py b/utils/enums.py index 22f91b8..a7982f8 100644 --- a/utils/enums.py +++ b/utils/enums.py @@ -50,9 +50,9 @@ class ParamTag(Enum): class PriorsCateg(Enum): - UNIFORM = 1 - GAUSSIAN = 2 - HALFGAUSS = 3 + UNIFORM = 1 + GAUSSIAN = 2 + LIMITEDGAUSS = 3 class MCMCSeedType(Enum): diff --git a/utils/fr.py b/utils/fr.py index 528557a..1d12fe5 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -97,9 +97,9 @@ def angles_to_fr(src_angles): spsi2 = SIN(psi)**2 cspi2 = 1. - spsi2 - x = sphi2*cspi2 - y = sphi2*spsi2 - z = cphi2 + x = float(abs(sphi2*cspi2)) + y = float(abs(sphi2*spsi2)) + z = float(abs(cphi2)) return x, y, z diff --git a/utils/gf.py b/utils/gf.py index 1998484..2c794d3 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -85,6 +85,11 @@ def steering_params(args): params.maxFitEnergy = 1E7 # GeV params.load_data_from_text_file = False + params.sampleToLoad = gf.sampleTag.MagicTau + params.use_legacy_selfveto_calculation = False + params.spline_hole_ice = False + params.spline_dom_efficiency = False + return params diff --git a/utils/llh.py b/utils/llh.py index f4404c4..0c1e97d 100644 --- a/utils/llh.py +++ b/utils/llh.py @@ -12,7 +12,8 @@ from __future__ import absolute_import, division from functools import partial import numpy as np -from scipy.stats import multivariate_normal, rv_continuous +import scipy +from scipy.stats import multivariate_normal, truncnorm from utils import fr as fr_utils from utils import gf as gf_utils @@ -20,10 +21,11 @@ from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg from utils.misc import enum_parse, gen_identifier, parse_bool -class Gaussian(rv_continuous): - """Gaussian for one dimension.""" - def _pdf(self, x, mu, sig): - return (1./np.sqrt(2*np.pi*sig**2))*np.exp(-((x-mu)**2)/(2*sig**2)) +def GaussianBoundedRV(loc=0., sigma=1., lower=-np.inf, upper=np.inf): + """Normalised gaussian bounded between lower and upper values""" + low, up = (lower - loc) / sigma, (upper - loc) / sigma + g = scipy.stats.truncnorm(loc=loc, scale=sigma, a=low, b=up) + return g def multi_gaussian(fr, fr_bf, sigma, offset=-320): @@ -69,13 +71,14 @@ def lnprior(theta, paramset): prior = 0 for param in paramset: if param.prior is PriorsCateg.GAUSSIAN: - prior += Gaussian().logpdf( - param.nominal_value, param.value, param.std - ) - elif param.prior is PriorsCateg.HALFGAUSS: - prior += Gaussian().logpdf( - param.nominal_value, param.value, param.std - ) + Gaussian().logcdf(1, param.value, param.std) + prior += GaussianBoundedRV( + loc=param.nominal_value, sigma=param.std + ).logpdf(param.value) + elif param.prior is PriorsCateg.LIMITEDGAUSS: + prior += GaussianBoundedRV( + loc=param.nominal_value, sigma=param.std, + lower=param.ranges[0], upper=param.ranges[1] + ).logpdf(param.value) return prior diff --git a/utils/plot.py b/utils/plot.py index 0529d5d..91b8b4e 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -15,7 +15,8 @@ from copy import deepcopy import numpy as np import numpy.ma as ma -from scipy import interpolate +from scipy.interpolate import splev, splprep +from scipy.ndimage.filters import gaussian_filter import matplotlib as mpl mpl.use('Agg') @@ -28,13 +29,18 @@ from matplotlib.lines import Line2D from matplotlib.patches import Patch from getdist import plots, mcsamples + import ternary +from ternary.heatmapping import polygon_generator + +import shapely.geometry as geometry from utils import misc as misc_utils from utils.enums import DataType, EnergyDependance from utils.enums import Likelihood, MixingScenario, ParamTag, StatCateg from utils.fr import angles_to_u, angles_to_fr + if os.path.isfile('./plot_llh/paper.mplstyle'): plt.style.use('./plot_llh/paper.mplstyle') elif os.environ.get('GOLEMSOURCEPATH') is not None: @@ -287,12 +293,13 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None): print 'data', data print 'data.shape', data.shape scales, statistic = ma.compress_rows(data).T - scales_rm = scales[1:] - statistic_rm = statistic[1:] - tck, u = interpolate.splprep([scales_rm, statistic_rm], s=0) - scales_rm, statistic_rm = interpolate.splev(np.linspace(0, 1, 1000), tck) - print 'scales_rm', scales_rm - print 'statistic_rm', statistic_rm + try: + tck, u = splprep([scales, statistic], s=0) + except: + return + sc, st = splev(np.linspace(0, 1, 10000), tck) + scales_rm = sc[sc >= scales[1]] + statistic_rm = st[sc >= scales[1]] min_idx = np.argmin(scales) null = statistic[min_idx] @@ -314,6 +321,10 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None): elif args.stat_method is StatCateg.FREQUENTIST: ax.set_ylabel(r'$-2\Delta {\rm LLH}$') + # ymin = np.round(np.min(reduced_ev) - 1.5) + # ymax = np.round(np.max(reduced_ev) + 1.5) + # ax.set_ylim((ymin, ymax)) + ax.plot(scales_rm, reduced_ev) ax.axhline(y=np.log(10**(3/2.)), color='red', alpha=1., linewidth=1.3) @@ -565,10 +576,14 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args): ax.axvline(x=ps, color='purple', alpha=1., linewidth=1.5) scales, statistic = ma.compress_rows(data[idim][isrc][ian]).T - scales_rm = scales[1:] - statistic_rm = statistic[1:] - tck, u = interpolate.splprep([scales_rm, statistic_rm], s=0) - scales_rm, statistic_rm = interpolate.splev(np.linspace(0, 1, 1000), tck) + try: + tck, u = splprep([scales, statistic], s=0) + except: + return + sc, st = splev(np.linspace(0, 1, 10000), tck) + scales_rm = sc[sc >= scales[1]] + statistic_rm = st[sc >= scales[1]] + min_idx = np.argmin(scales) null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: @@ -708,8 +723,8 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): for ian in xrange(len(data[idim][isrc])): print '=== an', ian scales, statistic = data[idim][isrc][ian].T - tck, u = interpolate.splprep([scales, statistic], s=0) - scales, statistic = interpolate.splev(np.linspace(0, 1, 1000), tck) + tck, u = splprep([scales, statistic], s=0) + scales, statistic = splev(np.linspace(0, 1, 1000), tck) min_idx = np.argmin(scales) null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: @@ -861,7 +876,7 @@ def plot_sens_corr_angle(data, outfile, outformat, args): print uni, c print len(uni) print np.unique(c) - + n = len(uni) assert len(np.unique(c)) == 1 c = c[0] @@ -872,7 +887,7 @@ def plot_sens_corr_angle(data, outfile, outformat, args): 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 @@ -891,7 +906,7 @@ def plot_sens_corr_angle(data, outfile, outformat, args): ax.set_xlim((mini, maxi)) ax.set_ylim(0, 1) ax.grid(b=False) - + x_v = reduced_pdat[0].round(decimals=4) y_v = reduced_pdat[1].round(decimals=4) uniques = np.unique(x_v) @@ -924,8 +939,8 @@ def plot_sens_corr_angle(data, outfile, outformat, args): 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) + tck, u = splprep([p_x, p_y], s=1e-5, per=True) + xi, yi = splev(np.linspace(0, 1, 1000), tck) except: xi, yi = p_x, p_y ax.fill(xi, yi, 'r', edgecolor='r', linewidth=1) @@ -953,6 +968,36 @@ def cmap_discretize(cmap, N): return mpl.colors.LinearSegmentedColormap(cmap.name + "_%d"%N, cdict, 1024) +def get_tax(ax, scale): + ax.set_aspect('equal') + + # Boundary and Gridlines + fig, tax = ternary.figure(ax=ax, scale=scale) + + # Draw Boundary and Gridlines + tax.boundary(linewidth=2.0) + tax.gridlines(color='grey', multiple=scale/5., linewidth=1.0, alpha=0.4, ls='--') + tax.gridlines(color='grey', multiple=scale/10., linewidth=0.5, alpha=0.4, ls=':') + + # Set Axis labels and Title + fontsize = 23 + tax.left_axis_label(r"$f_{\tau}$", fontsize=fontsize+8, offset=0.2, rotation=0) + tax.right_axis_label(r"$f_{\mu}$", fontsize=fontsize+8, offset=0.2, rotation=0) + tax.bottom_axis_label(r"$f_{e}$", fontsize=fontsize+8, position=(0.55, -0.20/2, 0.5), rotation=0) + + # Remove default Matplotlib axis + tax.get_axes().axis('off') + tax.clear_matplotlib_ticks() + + # Set ticks + ticks = np.linspace(0, 1, 6) + tax.ticks(ticks=ticks, locations=ticks*scale, axis='blr', linewidth=1, + offset=0.03, fontsize=fontsize, tick_formats='%.1f') + + tax._redraw_labels() + + return tax + def triangle_project(frs, llh, outfile, outformat, args, llh_paramset, fig_text): print 'Making triangle projection' fontsize = 23 @@ -976,8 +1021,8 @@ def triangle_project(frs, llh, outfile, outformat, args, llh_paramset, fig_text) mean = np.mean(llh) sig = np.std(llh) - min_scale = np.min(llh) - max_scale = np.max(mean+sig) + max_scale = np.max(llh) + min_scale = np.min(mean-sig) norm = mpl.colors.Normalize(vmin=min_scale, vmax=max_scale) colour = map(alp, map(cmap, map(norm, llh))) @@ -990,26 +1035,7 @@ def triangle_project(frs, llh, outfile, outformat, args, llh_paramset, fig_text) gs.update(hspace=0.3, wspace=0.15) ax = fig.add_subplot(gs[0]) - ax.set_aspect('equal') - - # Boundary and Gridlines - scale = 1 - fig, tax = ternary.figure(ax=ax, scale=scale) - - # Draw Boundary and Gridlines - tax.boundary(linewidth=2.0) - tax.gridlines(color='grey', multiple=scale/5., linewidth=1.0, alpha=0.4, ls='--') - tax.gridlines(color='grey', multiple=scale/10., linewidth=0.5, alpha=0.4, ls=':') - - # Set Axis labels and Title - fontsize = 23 - tax.left_axis_label(r"$f_{\tau}$", fontsize=fontsize+8, offset=0.2, rotation=0) - tax.right_axis_label(r"$f_{\mu}$", fontsize=fontsize+8, offset=0.2, rotation=0) - tax.bottom_axis_label(r"$f_{e}$", fontsize=fontsize+8, position=(0.55, -0.20/2, 0.5), rotation=0) - - # Remove default Matplotlib axis - tax.get_axes().axis('off') - tax.clear_matplotlib_ticks() + tax = get_tax(ax, scale=1) # Plot tax.scatter(frs, marker='o', s=0.1, color=colour) @@ -1034,13 +1060,92 @@ def triangle_project(frs, llh, outfile, outformat, args, llh_paramset, fig_text) cb.set_label(r'$LLH$', fontsize=fontsize+5, labelpad=20, horizontalalignment='center') - # Set ticks - tax.ticks(axis='blr', multiple=scale/5., linewidth=1, offset=0.03, - fontsize=fontsize, tick_formats='%.1f') - - tax._redraw_labels() - misc_utils.make_dir(outfile) for of in outformat: print 'Saving plot as {0}'.format(outfile+'.'+of) fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) + + +def heatmap(data, scale, vmin=None, vmax=None, style='triangular'): + for k, v in data.items(): + data[k] = np.array(v) + style = style.lower()[0] + if style not in ["t", "h", 'd']: + raise ValueError("Heatmap style must be 'triangular', 'dual-triangular', or 'hexagonal'") + + vertices_values = polygon_generator(data, scale, style) + + vertices = [] + for v, value in vertices_values: + vertices.append(v) + return vertices + + +def flavour_contour(frs, ax, nbins, coverage, linewidth=2, color='black'): + """Plot the flavour contour for a specified coverage.""" + # Histogram in flavour space + H, b = np.histogramdd( + (frs[:,0], frs[:,1], frs[:,2]), + bins=(nbins+1, nbins+1, nbins+1), range=((0, 1), (0, 1), (0, 1)) + ) + H = H / np.sum(H) + + # 3D smoothing + smoothing = 0.05 + H_s = gaussian_filter(H, sigma=smoothing) + + # Finding coverage + H_r = np.ravel(H_s) + H_rs = np.argsort(H_r)[::-1] + H_crs = np.cumsum(H_r[H_rs]) + thres = np.searchsorted(H_crs, coverage/100.) + mask_r = np.zeros(H_r.shape) + mask_r[H_rs[:thres]] = 1 + mask = mask_r.reshape(H_s.shape) + + # Get vertices inside covered region + binx = np.linspace(0, 1, nbins+1) + interp_dict = {} + for i in xrange(len(binx)): + for j in xrange(len(binx)): + for k in xrange(len(binx)): + if mask[i][j][k] == 1: + interp_dict[(i, j, k)] = H_s[i, j, k] + vertices = np.array(heatmap(interp_dict, nbins)) + points = vertices.reshape((len(vertices)*3, 2)) + + # Convex full to find points forming exterior bound + pc = geometry.MultiPoint(points) + polygon = pc.convex_hull + ex_cor = np.array(list(polygon.exterior.coords)) + + # Join points with a spline + tck, u = splprep([ex_cor.T[0], ex_cor.T[1]], s=0., per=1, k=1) + xi, yi = map(np.array, splev(np.linspace(0, 1, 300), tck)) + + # Spline again to smooth + tck, u = splprep([xi, yi], s=0.4, per=1, k=3) + xi, yi = map(np.array, splev(np.linspace(0, 1, 300), tck)) + ev_polygon = np.dstack((xi, yi))[0] + + def project_toflavour(p): + """Convert from cartesian to flavour space.""" + x, y = p + b = y / (np.sqrt(3)/2.) + a = x - b/2. + return [a, b, nbins-a-b] + + # Remove points interpolated outside flavour triangle + f_ev_polygon = np.array(map(project_toflavour, ev_polygon)) + xf, yf, zf = f_ev_polygon.T + mask = np.array((xf < 0) | (yf < 0) | (zf < 0) | (xf > nbins) | + (yf > nbins) | (zf > nbins)) + ev_polygon = np.dstack((xi[~mask], yi[~mask]))[0] + + # Plot + ax.plot( + ev_polygon.T[0], ev_polygon.T[1], linewidth=linewidth, color=color, + zorder=2, alpha=0.6, label=r'{0}\%'.format(int(coverage)) + ) + ax.scatter(points.T[0], points.T[1], marker='o', s=2, alpha=1, color=color, + zorder=3) |