From ca0ec62f2af59784b0ff2782037807b715b1a77b Mon Sep 17 00:00:00 2001 From: shivesh Date: Tue, 22 May 2018 08:16:43 -0500 Subject: Tue May 22 08:16:43 CDT 2018 --- plot_sens.py | 12 -------- sens.py | 21 +++++++------ submitter/mcmc_dag.py | 5 +-- submitter/sens_dag.py | 8 ++--- test/test_gf_simple.py | 17 ++++++----- utils/fr.py | 4 +-- utils/gf.py | 1 + utils/likelihood.py | 1 + utils/multinest.py | 3 +- utils/plot.py | 82 ++++++++++++++++++++++++++++++++++---------------- 10 files changed, 89 insertions(+), 65 deletions(-) diff --git a/plot_sens.py b/plot_sens.py index 3d1181a..6b82ebc 100755 --- a/plot_sens.py +++ b/plot_sens.py @@ -29,18 +29,6 @@ from utils.param import Param, ParamSet, get_paramsets from utils import multinest as mn_utils -import matplotlib as mpl -print mpl.rcParams.keys() -mpl.rcParams['text.usetex'] = True -# mpl.rcParams['text.latex.unicode'] = True -# mpl.rcParams['text.latex.preamble'] = r'\usepackage{cmbright}' -mpl.rcParams['font.family'] = 'sans-serif' -# mpl.rcParams['font.sans-serif'] = 'DejaVu Sans' -# mpl.rcParams['mathtext.fontset'] = 'stixsans' -# mpl.rcParams['mathtext.rm'] = 'DejaVu Sans' -# mpl.rcParams['mathtext.it'] = 'DejaVu Sans:italic' -# mpl.rcParams['mathtext.bf'] = 'DejaVu Sans:bold' - def define_nuisance(): """Define the nuisance parameters.""" diff --git a/sens.py b/sens.py index 58f8231..f8c8b56 100755 --- a/sens.py +++ b/sens.py @@ -273,17 +273,18 @@ def main(): print '|||| SCALE = {0:.0E}'.format(np.power(10, sc)) scale.value = sc if args.stat_method is StatCateg.BAYESIAN: - stat = mn_utils.mn_evidence( - mn_paramset = sens_paramset, - llh_paramset = llh_paramset, - asimov_paramset = asimov_paramset, - args = args, - fitter = fitter - ) - if stat is None: + try: + stat = mn_utils.mn_evidence( + mn_paramset = sens_paramset, + llh_paramset = llh_paramset, + asimov_paramset = asimov_paramset, + args = args, + fitter = fitter + ) + except: print 'Failed run, continuing' - raise - # continue + # raise + continue print '## Evidence = {0}'.format(stat) elif args.stat_method is StatCateg.FREQUENTIST: def fn(x): diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py index bad73ac..23586ac 100644 --- a/submitter/mcmc_dag.py +++ b/submitter/mcmc_dag.py @@ -25,7 +25,7 @@ GLOBAL_PARAMS = {} # MCMC GLOBAL_PARAMS.update(dict( - run_mcmc = 'True', + run_mcmc = 'False', burnin = 250, nsteps = 1000, nwalkers = 60, @@ -34,8 +34,9 @@ GLOBAL_PARAMS.update(dict( )) # FR -dimension = [3] +# dimension = [3, 6] # dimension = [4, 5, 7, 8] +dimension = [3, 4, 5, 6, 7, 8] GLOBAL_PARAMS.update(dict( threads = 1, binning = '6e4 1e7 20', diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py index 419ce2e..78034d2 100644 --- a/submitter/sens_dag.py +++ b/submitter/sens_dag.py @@ -27,9 +27,9 @@ GLOBAL_PARAMS = {} sens_eval_bin = 'true' # set to 'all' to run normally GLOBAL_PARAMS.update(dict( sens_run = 'True', - run_method = 'fixed_angle', # full, fixed_angle, corr_angle + run_method = 'corr_angle', # full, fixed_angle, corr_angle stat_method = 'bayesian', - sens_bins = 20, + sens_bins = 10, seed = 'None' )) @@ -43,10 +43,10 @@ GLOBAL_PARAMS.update(dict( # FR # dimension = [3] dimension = [3, 6] +# dimension = [4, 5, 7, 8] GLOBAL_PARAMS.update(dict( threads = 1, binning = '6e4 1e7 20', - # binning = '1e5 1e7 20', no_bsm = 'False', scale_region = "1E10", energy_dependance = 'spectral', @@ -65,7 +65,7 @@ GLOBAL_PARAMS.update(dict( # GolemFit GLOBAL_PARAMS.update(dict( ast = 'p2_0', - data = 'asimov' + data = 'real' )) # Plot diff --git a/test/test_gf_simple.py b/test/test_gf_simple.py index 3a6ebb5..1dcb392 100644 --- a/test/test_gf_simple.py +++ b/test/test_gf_simple.py @@ -5,28 +5,31 @@ import matplotlib.pyplot as plt import GolemFitPy as gf -FASTMODE = True +FASTMODE = False dp = gf.DataPaths() npp = gf.NewPhysicsParams() -sp = gf.SteeringParams(gf.sampleTag.HESE) +sp = gf.SteeringParams(gf.sampleTag.MagicTau) sp.quiet = False if FASTMODE: sp.fastmode = True sp.frequentist = True +sp.load_data_from_text_file = False golem = gf.GolemFit(dp, sp, npp) -fp = gf.FitParameters(gf.sampleTag.HESE) +fp = gf.FitParameters(gf.sampleTag.MagicTau) fp.astroFlavorAngle1 = 4./9. -fp.astroFlavorAngle2 = 0 +fp.astroFlavorAngle2 = 0. golem.SetupAsimov(fp) -fp_sh = gf.FitParameters(gf.sampleTag.HESE) -fp_sh.astroFlavorAngle1 = 0.36 -fp_sh.astroFlavorAngle2 = -0.57 +fp_sh = gf.FitParameters(gf.sampleTag.MagicTau) +# fp_sh.astroFlavorAngle1 = 0.36 +# fp_sh.astroFlavorAngle2 = -0.57 +fp_sh.astroFlavorAngle1 = 1. +fp_sh.astroFlavorAngle2 = 0. print 'Eval fp = {0}'.format(golem.EvalLLH(fp)) diff --git a/utils/fr.py b/utils/fr.py index 51864a3..cb830ea 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -144,9 +144,9 @@ def angles_to_u(bsm_angles): c23 = COS(t23) s23 = SIN(t23) - p1 = np.array([[1 , 0 , 0] , [0 , c23 , s23] , [0 , -s23 , c23]] , dtype=CDTYPE) + p1 = np.array([[1 , 0 , 0] , [0 , c23 , s23] , [0 , -s23 , c23]] , dtype=CDTYPE) p2 = np.array([[c13 , 0 , s13*EXP(-1j*dcp)] , [0 , 1 , 0] , [-s13*EXP(1j*dcp) , 0 , c13]] , dtype=CDTYPE) - p3 = np.array([[c12 , s12 , 0] , [-s12 , c12 , 0] , [0 , 0 , 1]] , dtype=CDTYPE) + p3 = np.array([[c12 , s12 , 0] , [-s12 , c12 , 0] , [0 , 0 , 1]] , dtype=CDTYPE) u = np.dot(np.dot(p1, p2), p3) return u diff --git a/utils/gf.py b/utils/gf.py index 79f1dd0..06a6125 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -92,6 +92,7 @@ def setup_fitter(args, asimov_paramset): def get_llh(fitter, params): + # print 'params', params fitparams = gf.FitParameters(gf.sampleTag.MagicTau) for parm in params: fitparams.__setattr__(parm.name, float(parm.value)) diff --git a/utils/likelihood.py b/utils/likelihood.py index ed598c8..ee4ec6e 100644 --- a/utils/likelihood.py +++ b/utils/likelihood.py @@ -164,6 +164,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): fr = averaged_measured_flux / np.sum(averaged_measured_flux) flavour_angles = fr_utils.fr_to_angles(fr) + # print 'flavour_angles', map(float, flavour_angles) for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): param.value = flavour_angles[idx] diff --git a/utils/multinest.py b/utils/multinest.py index 387ff90..63437bc 100644 --- a/utils/multinest.py +++ b/utils/multinest.py @@ -84,7 +84,7 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter): ) make_dir(prefix) print 'Running evidence calculation for {0}'.format(prefix) - result = run( + run( LogLikelihood = lnProbEval, Prior = CubePrior, n_dims = n_params, @@ -95,7 +95,6 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter): resume = False, verbose = True ) - if result is None: return None analyser = analyse.Analyzer( outputfiles_basename=prefix, n_params=n_params diff --git a/utils/plot.py b/utils/plot.py index b63126b..4dd8cc4 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -13,6 +13,8 @@ import os from copy import deepcopy import numpy as np +from scipy import interpolate + import matplotlib as mpl mpl.use('Agg') from matplotlib import rc @@ -27,7 +29,7 @@ from utils.enums import Likelihood, ParamTag, StatCateg from utils.fr import angles_to_u, angles_to_fr rc('text', usetex=False) -rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) +rc('font', **{'family':'serif', 'serif':['DejaVu Sans'], 'size':18}) def centers(x): @@ -117,17 +119,18 @@ def gen_figtext(args): 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}]\nIC ' \ - 'observed flavour ratio = [{3}, {4}, ' \ - '{5}]\nDimension = {6}'.format( - sr1, sr2, sr3, mr1, mr2, mr3, args.dimension, - int(args.energy) - ) + t += 'Source flavour ratio = [{0}, {1}, {2}]'.format(sr1, sr2, sr3) + if args.data is DataType.ASIMOV: + t += '\nIC observed flavour ratio = [{0}, {1}, {2}]'.format( + mr1, mr2, mr3 + ) + t += '\nDimension = {0}'.format(args.dimension) else: - t += 'IC observed flavour ratio = [{0}, {1}, ' \ - '{2}]\nDimension = {3}'.format( - mr1, mr2, mr3, args.dimension, int(args.energy) - ) + if args.data is DataType.ASIMOV: + 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: @@ -231,7 +234,7 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset): g.export(outfile+'_elements.'+of) -def myround(x, base=1, up=False, down=False): +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)) @@ -241,12 +244,17 @@ def myround(x, base=1, up=False, down=False): def plot_statistic(data, outfile, outformat, args, scale_param, label=None): """Make MultiNest factor or LLH value plot.""" print 'Making Statistic plot' - fig_text = gen_figtext(args) + fig_text = '\n' + gen_figtext(args) if label is not None: fig_text += '\n' + label print 'data', data print 'data.shape', data.shape scales, statistic = data.T + tck, u = interpolate.splprep([scales, statistic], s=0) + scales, statistic = interpolate.splev(np.linspace(0, 1, 1000), tck) + print 'scales', scales + print 'statistic', statistic + min_idx = np.argmin(scales) null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: @@ -258,21 +266,31 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None): ax = fig.add_subplot(111) ax.set_xlim(np.log10(args.scale_region)) - ax.set_xlabel('$'+scale_param.tex+'$') + ax.set_xlabel(r'${\mathrm {log}}_{10} \left (\Lambda^{-1}' + \ + get_units(args.dimension) +r'\right )$', fontsize=16) if args.stat_method is StatCateg.BAYESIAN: - ax.set_ylabel(r'Bayes Factor') + ax.set_ylabel(r'log(Bayes Factor)') elif args.stat_method is StatCateg.FREQUENTIST: ax.set_ylabel(r'$-2\Delta {\rm LLH}$') ax.plot(scales, reduced_ev) + ax.axhline(y=np.log(10**(3/2.)), 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) for xmaj in ax.xaxis.get_majorticklocs(): ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.3, linewidth=1) + if args.data is DataType.REAL: + fig.text(0.8, 0.14, 'IceCube Preliminary', color='red', fontsize=9, + ha='center', va='center') + elif args.data is DataType.ASIMOV: + fig.text(0.8, 0.14, 'IceCube Simulation', color='red', fontsize=9, + ha='center', va='center') + at = AnchoredText( - fig_text, prop=dict(size=10), frameon=True, loc=2 + fig_text, prop=dict(size=10), frameon=True, loc=4 ) at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5") ax.add_artist(at) @@ -374,9 +392,13 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): fig = plt.figure(figsize=(7, 5)) ax = fig.add_subplot(111) ax.set_xlim(0, len(xticks)+1) - ax.set_xticklabels([''] + xticks + ['']) - ax.set_xlabel(r'BSM operator angle') - ax.set_ylabel(r'${\rm log}_{10} \left (\Lambda^{-1}' + get_units(dim) +r'\right )$') + ax.set_xticklabels([''] + xticks + [''], fontsize=14) + ax.set_xlabel(r'BSM operator angle', fontsize=16) + ax.set_ylabel(r'${\mathrm {log}}_{10} \left (\Lambda^{-1}' + \ + get_units(dim) +r'\right )$', fontsize=16) + + # ax.tick_params(axis='x', labelsize=11) + ax.tick_params(axis='y', labelsize=14) for isrc in xrange(len(data[idim])): src = args.source_ratios[isrc] @@ -386,6 +408,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) min_idx = np.argmin(scales) null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: @@ -400,24 +424,30 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): dim, src, reduced_ev ) continue + if reduced_ev[-1] < np.log(10**(3/2.)) - 0.1: + print 'Peaked contour does not exclude large scales! For ' \ + 'DIM {0} FRS{1} reduced_ev {2}!'.format( + dim, src, reduced_ev + ) + continue arr_len = dim-2 lim = al[0] print 'limit = {0}'.format(lim) label = '{0} : {1} : {2}'.format(*misc_utils.solve_ratio(src)) - if lim < yranges[0]: yranges[0] = lim + if lim < yranges[0]: yranges[0] = lim-arr_len if lim > yranges[1]: yranges[1] = lim+arr_len+2 # if lim > yranges[1]: yranges[1] = lim xoff = 0.15 line = plt.Line2D( - (ian+1-xoff, ian+1+xoff), (lim, lim), lw=2., color=colour[isrc], label=label + (ian+1-xoff, ian+1+xoff), (lim, lim), lw=2.5, color=colour[isrc], label=label ) ax.add_line(line) if len(legend_handles) < isrc+1: legend_handles.append(line) x_offset = isrc*xoff/2. - xoff/2. ax.annotate( - s='', xy=(ian+1+x_offset, lim-0.01), xytext=(ian+1+x_offset, lim+arr_len), - arrowprops={'arrowstyle': '<-', 'lw': 2., 'color':colour[isrc]} + s='', xy=(ian+1+x_offset, lim-0.02), xytext=(ian+1+x_offset, lim+arr_len), + arrowprops={'arrowstyle': '<-', 'lw': 1.5, 'color':colour[isrc]} ) try: @@ -438,14 +468,14 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5") ax.add_artist(at) - fig.text(0.42, 0.8, 'Excluded', color='red', fontsize=20, ha='center', + fig.text(0.42, 0.8, r'Excluded', color='red', fontsize=16, ha='center', va='center', fontweight='bold') if args.data is DataType.REAL: - fig.text(0.805, 0.14, 'IceCube Preliminary', color='red', fontsize=11, + fig.text(0.805, 0.14, 'IceCube Preliminary', color='red', fontsize=9, ha='center', va='center') elif args.data is DataType.ASIMOV: - fig.text(0.805, 0.14, 'IceCube Simulation', color='red', fontsize=11, + fig.text(0.805, 0.14, 'IceCube Simulation', color='red', fontsize=9, ha='center', va='center') for ymaj in ax.yaxis.get_majorticklocs(): -- cgit v1.2.3