diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rwxr-xr-x | fr.py | 93 | ||||
| -rw-r--r-- | mnrun/plot.py | 23 | ||||
| -rw-r--r-- | mnrun/test.png | bin | 0 -> 64972 bytes | |||
| -rwxr-xr-x | sens.py | 146 | ||||
| -rw-r--r-- | submitter/make_dag.py | 26 | ||||
| -rw-r--r-- | utils/enums.py | 1 | ||||
| -rw-r--r-- | utils/fr.py | 17 | ||||
| -rw-r--r-- | utils/gf.py | 11 | ||||
| -rw-r--r-- | utils/likelihood.py | 10 | ||||
| -rw-r--r-- | utils/plot.py | 6 |
11 files changed, 306 insertions, 28 deletions
@@ -6,3 +6,4 @@ plots/ *.pdf dagman_FR.submit* submitter/logs/ +mnrun_* @@ -22,7 +22,7 @@ from utils import likelihood as llh_utils from utils import mcmc as mcmc_utils from utils import misc as misc_utils from utils.enums import EnergyDependance, Likelihood, MCMCSeedType, ParamTag -from utils.fr import MASS_EIGENVALUES, normalise_fr +from utils.fr import MASS_EIGENVALUES, normalise_fr, fr_to_angles from utils.misc import enum_parse, Param, ParamSet from utils.plot import plot_argparse, chainer_plot @@ -64,10 +64,10 @@ def get_paramsets(args): 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='astroENorm' , value=args.measured_ratio[0], ranges=[0., 1.], std=0.2, tag=tag), - Param(name='astroMuNorm' , value=args.measured_ratio[1], ranges=[0., 1.], std=0.2, tag=tag), - Param(name='astroTauNorm', value=args.measured_ratio[2], ranges=[0., 1.], std=0.2, tag=tag) + 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) @@ -110,6 +110,10 @@ def process_args(args): raise NotImplementedError( '--fix-mixing and --fix-mixing-almost cannot be used together' ) + if args.bayes_factor and args.fix_scale: + raise NotImplementedError( + '--bayes-factor and --fix-scale cannot be used together' + ) args.measured_ratio = normalise_fr(args.measured_ratio) if args.fix_source_ratio: @@ -164,10 +168,22 @@ def parse_args(): help='Spectral index for spectral energy dependance' ) parser.add_argument( - '--binning', default=[1e4, 1e7, 10], type=float, nargs=3, + '--binning', default=[1e4, 1e7, 5], type=float, nargs=3, help='Binning for spectral energy dependance' ) parser.add_argument( + '--bayes-factor', type=misc_utils.parse_bool, default='False', + help='Make the bayes factor plot for the scale' + ) + parser.add_argument( + '--bayes-bins', type=int, default=100, + help='Number of bins for the Bayes factor plot' + ) + parser.add_argument( + '--bayes-output', type=str, default='./mnrun/', + help='Folder to store MultiNest evaluations' + ) + parser.add_argument( '--source-ratio', type=int, nargs=3, default=[2, 1, 0], help='Set the source flavour ratio for the case when you want to fix it' ) @@ -239,7 +255,6 @@ def main(): datapaths = gf.DataPaths() sparams = gf_utils.steering_params(args) npp = gf.NewPhysicsParams() - fitter = gf.GolemFit(datapaths, sparams, npp) gf_utils.set_up_as(fitter, asimov_paramset) # fitter.WriteCompact() @@ -274,7 +289,6 @@ def main(): ) mcmc_utils.save_chains(samples, outfile) - print "Making triangle plots" chainer_plot( infile = outfile+'.npy', outfile = outfile[:5]+outfile[5:].replace('data', 'plots'), @@ -282,6 +296,71 @@ def main(): args = args, mcmc_paramset = mcmc_paramset ) + + if args.bayes_factor: + # TODO(shivesh) + import pymultinest + from pymultinest.solve import solve + from pymultinest.watch import ProgressPrinter + + if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT: + datapaths = gf.DataPaths() + sparams = gf_utils.steering_params(args) + npp = gf.NewPhysicsParams() + fitter = gf.GolemFit(datapaths, sparams, npp) + gf_utils.set_up_as(fitter, asimov_paramset) + else: + fitter = None + + sc_range = mcmc_paramset.from_tag(ParamTag.SCALE)[0].ranges + scales = np.linspace( + sc_range[0], sc_range[1], args.bayes_bins+1 + ) + + p = mcmc_paramset.from_tag(ParamTag.SCALE, invert=True) + prior_ranges = p.ranges + n_params = len(p) + hypo_paramset = asimov_paramset + + def CubePrior(cube, ndim, nparams): + # default are uniform priors + return ; + + arr = [] + for s in scales: + theta = np.zeros(n_params) + def lnProb(cube, ndim, nparams): + for i in range(ndim): + prange = prior_ranges[i][1] - prior_ranges[i][0] + theta[i] = prange*cube[i] + prior_ranges[i][0] + theta_ = np.concatenate([theta, [s]]) + return llh_utils.triangle_llh( + theta=theta_, + args=args, + asimov_paramset=asimov_paramset, + mcmc_paramset=mcmc_paramset, + fitter=fitter + ) + + prefix = 'mnrun_{0:.0E}'.format(np.power(10, s)) + print 'begin running evidence calculation for {0}'.format(prefix) + result = pymultinest.run( + LogLikelihood=lnProb, Prior=CubePrior, n_dims=n_params, + outputfiles_basename=prefix#, verbose=True + ) + print 'end running evidence calculation for {0}'.format(prefix) + + print 'begin analysis for {0}'.format(prefix) + analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params) + a_lnZ = analyzer.get_stats()['global evidence'] + print 'end analysis for {0}'.format(prefix) + print 'Evidence = {0}'.format(a_lnZ) + + arr.append([s, a_lnZ]) + out = args.bayes_output+'fr_evidence_'+outfile[2:]+'.npy' + misc_utils.make_dir(out) + np.save(out, np.array(arr)) + print "DONE!" diff --git a/mnrun/plot.py b/mnrun/plot.py new file mode 100644 index 0000000..e010257 --- /dev/null +++ b/mnrun/plot.py @@ -0,0 +1,23 @@ +import numpy as np +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import pyplot as plt +from matplotlib import rc + +rc('text', usetex=False) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +arr = np.load('fr_evidence_test_050_050_000_0100_sfr_033_066_000_DIM3_single_scale.npy') + +print arr +ev = arr.T[1] +null = ev[0] +print null +re_z = -(ev - null) +print re_z + +plt.plot(arr.T[0], re_z) +plt.xlabel(r'${\rm log}_{10} \Lambda$') +plt.ylabel(r'Bayes Factor') +# plt.axhline(np.power(10, 1.5), color='grey', linestyle='-', alpha=0.4) +plt.savefig('./test.png', bbox_inches='tight', dpi=150) diff --git a/mnrun/test.png b/mnrun/test.png Binary files differnew file mode 100644 index 0000000..3f68e17 --- /dev/null +++ b/mnrun/test.png @@ -0,0 +1,146 @@ +#! /usr/bin/env python +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : April 10, 2018 + +""" +HESE BSM flavour ratio analysis script +""" + +from __future__ import absolute_import, division + +import numpy as np +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import pyplot as plt +from matplotlib import rc + +import GolemFitPy as gf + +import fr +from utils import gf as gf_utils +from utils import likelihood as llh_utils +from utils import misc as misc_utils +from utils.enums import Likelihood, ParamTag + +rc('text', usetex=False) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + + +RUN = True + + +z = 0+1E-6 +scenarios = [ + [np.sin(np.pi/2.)**2, z, z, z], + [z, np.cos(np.pi/2.)**4, z, z], + [z, z, np.sin(np.pi/2.)**2, z], +] +xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$'] + +def fit_flags(flag_dict): + flags = gf.FitParametersFlag() + for key in flag_dict.iterkeys(): + flags.__setattr__(key, flag_dict[key]) + return flags + +default_flags = { + # False means it's not fixed in minimization + 'astroFlavorAngle1' : True, + 'astroFlavorAngle2' : True, + # 'astroENorm' : True, + # 'astroMuNorm' : True, + # 'astroTauNorm' : True, + 'convNorm' : False, + 'promptNorm' : False, + 'muonNorm' : False, + 'astroNorm' : False, + 'astroParticleBalance' : True, + 'astroDeltaGamma' : False, + 'cutoffEnergy' : True, + 'CRDeltaGamma' : False, + 'piKRatio' : False, + 'NeutrinoAntineutrinoRatio' : True, + 'darkNorm' : True, + 'domEfficiency' : True, + 'holeiceForward' : True, + 'anisotropyScale' : True, + 'astroNormSec' : True, + 'astroDeltaGammaSec' : True +} + + +def main(): + args = fr.parse_args() + args.likelihood = Likelihood.GF_FREQ + fr.process_args(args) + misc_utils.print_args(args) + + asimov_paramset, mcmc_paramset = fr.get_paramsets(args) + outfile = misc_utils.gen_outfile_name(args) + print '== {0:<25} = {1}'.format('outfile', outfile) + + asimov_paramset = asimov_paramset.from_tag(ParamTag.BESTFIT) + mcmc_paramset = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True) + + sc_range = mcmc_paramset.from_tag(ParamTag.SCALE)[0].ranges + scan_scales = np.linspace(sc_range[0], sc_range[1], 100) + print 'scan_scales', scan_scales + + if RUN: + datapaths = gf.DataPaths() + sparams = gf_utils.steering_params(args) + npp = gf.NewPhysicsParams() + fitter = gf.GolemFit(datapaths, sparams, npp) + fitter.SetFitParametersFlag(fit_flags(default_flags)) + gf_utils.set_up_as(fitter, asimov_paramset) + + x = [] + y = [] + mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES) + for idx, scen in enumerate(scenarios): + scales = [] + llhs = [] + for yidx, an in enumerate(mm_angles): + an.value = scen[yidx] + for sc in scan_scales: + theta = scen + [sc] + print 'theta', theta + llh = llh_utils.triangle_llh( + theta=theta, args=args, asimov_paramset=asimov_paramset, + mcmc_paramset=mcmc_paramset, fitter=fitter + ) + print 'llh', llh + scales.append(sc) + llhs.append(llh) + min_llh = np.min(llhs) + delta_llh = 2*(np.array(llhs) - min_llh) + print 'scales', scales + print 'delta_llh', delta_llh + + n_arr = [] + for i, d in enumerate(delta_llh): + # 90% for 1 DOF + if d < 2.71: + x.append(idx) + y.append(scales[i]) + + np.save('sens.npy', np.array([x, y])) + else: + x, y = np.load('sens.npy') + + plt.plot(x, y, linewidth=0., marker='.') + plt.xticks(range(len(scenarios)), xticks) + plt.xlim(-1, len(scenarios)) + plt.ylim(sc_range[0], sc_range[1]) + plt.ylabel(r'${\rm log}_{10}\Lambda / GeV$') + plt.savefig('./sens.png', bbox_inches='tight', dpi=150) + + +main.__doc__ = __doc__ + + +if __name__ == '__main__': + main() + diff --git a/submitter/make_dag.py b/submitter/make_dag.py index 15d195a..735d213 100644 --- a/submitter/make_dag.py +++ b/submitter/make_dag.py @@ -16,17 +16,17 @@ full_scan_mfr = [ ] fix_sfr_mfr = [ - # (1, 1, 1, 1, 0, 0), - # (1, 1, 1, 0, 1, 0), - # (1, 1, 1, 0, 0, 1), + (1, 1, 1, 1, 0, 0), + (1, 1, 1, 0, 1, 0), + (1, 1, 1, 0, 0, 1), (1, 1, 1, 1, 2, 0), (1, 1, 0, 0, 1, 0), - # (1, 1, 0, 1, 2, 0), - # (1, 1, 0, 1, 0, 0), - # (1, 0, 0, 1, 0, 0), - # (0, 1, 0, 0, 1, 0), - # (1, 2, 0, 0, 1, 0), - # (1, 2, 0, 1, 2, 0) + (1, 1, 0, 1, 2, 0), + (1, 1, 0, 1, 0, 0), + (1, 0, 0, 1, 0, 0), + (0, 1, 0, 0, 1, 0), + (1, 2, 0, 0, 1, 0), + (1, 2, 0, 1, 2, 0) ] # MCMC @@ -39,7 +39,7 @@ threads = 4 mcmc_seed_type = 'uniform' # FR -dimension = [3, 6] +dimension = [4, 5, 7, 8] energy = [1e6] likelihood = 'golemfit' no_bsm = 'False' @@ -48,9 +48,9 @@ scale = "1E-20 1E-30" scale_region = "1E10" energy_dependance = 'spectral' spectral_index = -2 -binning = [1e4, 1e7, 10] +binning = [1e4, 1e7, 5] fix_mixing = 'False' -fix_mixing_almost = 'True' +fix_mixing_almost = 'False' # Likelihood likelihood = 'golemfit' @@ -70,7 +70,7 @@ data = 'real' plot_angles = 'True' plot_elements = 'False' -outfile = 'dagman_FR_fix_mixing_almost.submit' +outfile = 'dagman_FR.submit' golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit' condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/submit.sub' diff --git a/utils/enums.py b/utils/enums.py index 4637eeb..90ca17d 100644 --- a/utils/enums.py +++ b/utils/enums.py @@ -30,6 +30,7 @@ class Likelihood(Enum): FLAT = 1 GAUSSIAN = 2 GOLEMFIT = 3 + GF_FREQ = 4 class ParamTag(Enum): diff --git a/utils/fr.py b/utils/fr.py index 7f9d855..8f71f78 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -194,6 +194,23 @@ def normalise_fr(fr): return np.array(fr) / float(np.sum(fr)) +def fr_to_angles(ratios): + """Convert from flavour ratio into the angular projection of the flavour + ratios. + + Parameters + ---------- + TODO(shivesh) + """ + fr0, fr1, fr2 = normalise_fr(ratios) + sphi4 = (fr2 - 1.0)**2 + if (fr2 - 1.0) == 0: + c2psi = 0 + else: + c2psi = (fr1*2.0 + fr2 - 1.0) * (fr2 - 1.0) + return sphi4, c2psi + + NUFIT_U = angles_to_u((0.307, (1-0.02195)**2, 0.565, 3.97935)) """NuFIT mixing matrix (s_12^2, c_13^4, s_23^2, dcp)""" diff --git a/utils/gf.py b/utils/gf.py index ebc8538..e575094 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -44,14 +44,21 @@ def set_up_as(fitter, params): def get_llh(fitter, params): fitparams = gf.FitParameters(gf.sampleTag.HESE) - # print params for parm in params: fitparams.__setattr__(parm.name, parm.value) llh = -fitter.EvalLLH(fitparams) - # print '=== llh = {0}'.format(llh) return llh +def get_llh_freq(fitter, params): + print 'setting to {0}'.format(params) + fitparams = gf.FitParameters(gf.sampleTag.HESE) + for parm in params: + fitparams.__setattr__(parm.name, parm.value) + fitter.SetFitParametersSeed(fitparams) + return -fitter.MinLLH().likelihood + + def data_distributions(fitter): hdat = fitter.GetDataDistribution() binedges = np.asarray([fitter.GetZenithBinsData(), fitter.GetEnergyBinsData()]) diff --git a/utils/likelihood.py b/utils/likelihood.py index fe4301d..c1208ab 100644 --- a/utils/likelihood.py +++ b/utils/likelihood.py @@ -101,7 +101,7 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): elif args.energy_dependance is EnergyDependance.SPECTRAL: mf_perbin = [] for i_sf, sf_perbin in enumerate(source_flux): - u = fr_utils.params_to_BSMu( + u = fr_utils.params_to_BSMu( theta = bsm_angles, dim = args.dimension, energy = args.energy, @@ -119,10 +119,9 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): intergrated_measured_flux fr = averaged_measured_flux / np.sum(averaged_measured_flux) + flavour_angles = fr_utils.fr_to_angles(fr) for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): - param.value = fr[idx] - - # print 'hypo_paramset', hypo_paramset + param.value = flavour_angles[idx] if args.likelihood is Likelihood.FLAT: return 1. @@ -131,6 +130,9 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): return gaussian_llh(fr, fr_bf, args.sigma_ratio) elif args.likelihood is Likelihood.GOLEMFIT: return gf_utils.get_llh(fitter, hypo_paramset) + elif args.likelihood is Likelihood.GF_FREQ: + 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) diff --git a/utils/plot.py b/utils/plot.py index c70ffec..4c5ff1c 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -171,6 +171,7 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): ranges = mcmc_paramset.ranges if args.plot_angles: + print "Making triangle plots" Tchain = raw g = plot_Tchain(Tchain, axes_labels, ranges) @@ -182,8 +183,8 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): itv = interval(Tchain[:,i_ax_1], percentile=90.) for l in itv: ax_2.axvline(l, color='gray', ls='--') - ax_2.set_title(r'${0:.2f}_{{-{1:.2f}}}^{{+{2:.2f}}}$'.format( - itv[1], itv[0], itv[2] + ax_2.set_title(r'${0:.2f}_{{{1:.2f}}}^{{+{2:.2f}}}$'.format( + itv[1], itv[0]-itv[1], itv[2]-itv[1] ), fontsize=10) # if not args.fix_mixing: @@ -198,6 +199,7 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): g.export(outfile+'_angles.'+of) if args.plot_elements: + print "Making triangle plots" if args.fix_mixing_almost: raise NotImplementedError nu_index = mcmc_paramset.from_tag(ParamTag.NUISANCE, index=True) |
