aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xfr.py5
-rwxr-xr-xfr_edit.py581
-rwxr-xr-xsens.py74
-rw-r--r--submitter/mcmc_dag.py24
-rw-r--r--submitter/sens_dag.py31
-rw-r--r--submitter/sens_submit.sub2
-rw-r--r--utils/fr.py3
-rw-r--r--utils/gf.py3
-rw-r--r--utils/likelihood.py6
-rw-r--r--utils/misc.py1
-rw-r--r--utils/multinest.py2
-rw-r--r--utils/plot.py5
12 files changed, 87 insertions, 650 deletions
diff --git a/fr.py b/fr.py
index 7558c8e..971ee6e 100755
--- a/fr.py
+++ b/fr.py
@@ -81,6 +81,10 @@ def process_args(args):
args.binning = np.logspace(
np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1
)
+ if args.likelihood is Likelihood.GOLEMFIT:
+ print 'GolemFit selected with spectral index energy dependance, ' \
+ 'will attempt to use the astroDeltaGamma systematic to fold ' \
+ 'in the spectral index.'
if not args.fix_scale:
args.scale = fr_utils.estimate_scale(args)
@@ -174,4 +178,3 @@ main.__doc__ = __doc__
if __name__ == '__main__':
main()
-
diff --git a/fr_edit.py b/fr_edit.py
deleted file mode 100755
index df3b43b..0000000
--- a/fr_edit.py
+++ /dev/null
@@ -1,581 +0,0 @@
-#! /usr/bin/env python
-# author : S. Mandalia
-# s.p.mandalia@qmul.ac.uk
-#
-# date : March 17, 2018
-
-"""
-HESE BSM flavour ratio analysis script
-"""
-
-from __future__ import absolute_import, division
-
-import os
-import argparse
-from functools import partial
-
-import numpy as np
-
-import GolemFitPy as gf
-
-from utils import gf as gf_utils
-from utils import likelihood as llh_utils
-from utils import mcmc as mcmc_utils
-from utils import misc as misc_utils
-from utils import plot as plot_utils
-from utils.enums import EnergyDependance, Likelihood, MCMCSeedType, ParamTag
-from utils.fr import MASS_EIGENVALUES, normalise_fr, fr_to_angles
-from utils.misc import enum_parse, Param, ParamSet
-
-
-def define_nuisance():
- """Define the nuisance parameters to scan over with default values,
- ranges and sigma.
- """
- tag = ParamTag.NUISANCE
- nuisance = ParamSet(
- 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)
- )
- return nuisance
-
-
-def nuisance_argparse(parser):
- nuisance_paramset = define_nuisance()
- for parm in nuisance_paramset:
- parser.add_argument(
- '--'+parm.name, type=float, default=parm.value,
- help=parm.name+' to inject'
- )
-
-
-def get_paramsets(args):
- """Make the paramsets for generating the Asmimov MC sample and also running
- the MCMC.
- """
- asimov_paramset = []
- mcmc_paramset = []
- if args.likelihood == Likelihood.GOLEMFIT:
- nuisance_paramset = define_nuisance()
- asimov_paramset.extend(nuisance_paramset.params)
- mcmc_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
- mcmc_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., 2*np.pi], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag)
- ])
- if args.fix_mixing_almost:
- tag = ParamTag.MMANGLES
- mcmc_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
- mcmc_paramset.append(
- Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
- tex=r'{\rm log}_{10}\Lambda'+plot_utils.get_units(args.dimension), tag=tag)
- )
- if not args.fix_source_ratio:
- tag = ParamTag.SRCANGLES
- mcmc_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)
- ])
- mcmc_paramset = ParamSet(mcmc_paramset)
- return asimov_paramset, mcmc_paramset
-
-
-def process_args(args):
- """Process the input args."""
- if args.fix_mixing and args.fix_scale:
- raise NotImplementedError('Fixed mixing and scale not implemented')
- if args.fix_mixing and args.fix_mixing_almost:
- raise NotImplementedError(
- '--fix-mixing and --fix-mixing-almost cannot be used together'
- )
- if args.run_bayes_factor and args.fix_scale:
- raise NotImplementedError(
- '--run-bayes-factor and --fix-scale cannot be used together'
- )
-
- args.measured_ratio = normalise_fr(args.measured_ratio)
- if args.fix_source_ratio:
- args.source_ratio = normalise_fr(args.source_ratio)
-
- if args.energy_dependance is EnergyDependance.SPECTRAL:
- args.binning = np.logspace(
- np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1
- )
-
- if not args.fix_scale:
- if args.energy_dependance is EnergyDependance.MONO:
- args.scale = np.power(
- 10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \
- np.log10(args.energy**(args.dimension-3)) + 6
- )
- elif args.energy_dependance is EnergyDependance.SPECTRAL:
- args.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))
- )
- )
- """Estimate the scale of when to expect the BSM term to contribute."""
- args.scale_region = (args.scale/args.scale_region, args.scale*args.scale_region)
-
- if args.bayes_eval_bin.lower() == 'all':
- args.bayes_eval_bin = None
- else:
- args.bayes_eval_bin = map(int, args.bayes_eval_bin.split())
-
-
-def parse_args(args=None):
- """Parse command line arguments"""
- parser = argparse.ArgumentParser(
- description="BSM flavour ratio analysis",
- formatter_class=misc_utils.SortingHelpFormatter,
- )
- 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(
- '--sigma-ratio', type=float, default=0.01,
- help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH'
- )
- parser.add_argument(
- '--fix-source-ratio', type=misc_utils.parse_bool, default='False',
- help='Fix the source flavour ratio'
- )
- 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(
- '--run-bayes-factor', type=misc_utils.parse_bool, default='False',
- help='Make the bayes factor plot for the scale'
- )
- parser.add_argument(
- '--run-angles-limit', type=misc_utils.parse_bool, default='False',
- help='Make the limit vs BSM angles plot'
- )
- parser.add_argument(
- '--run-angles-correlation', type=misc_utils.parse_bool, default='False',
- help='Make the limit vs BSM angles plot'
- )
- parser.add_argument(
- '--bayes-bins', type=int, default=10,
- help='Number of bins for the Bayes factor plot'
- )
- parser.add_argument(
- '--bayes-eval-bin', type=str, default='all',
- help='Which bin to evalaute for Bayes factor plot'
- )
- parser.add_argument(
- '--bayes-live-points', type=int, default=400,
- help='Number of live points for MultiNest evaluations'
- )
- parser.add_argument(
- '--bayes-tolerance', type=float, default=0.5,
- help='Tolerance for MultiNest'
- )
- parser.add_argument(
- '--bayes-output', type=str, default='./mnrun/',
- help='Folder to store MultiNest evaluations'
- )
- parser.add_argument(
- '--angles-lim-output', type=str, default='./mnrun/',
- help='Folder to store MultiNest evaluations'
- )
- parser.add_argument(
- '--angles-corr-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'
- )
- parser.add_argument(
- '--no-bsm', type=misc_utils.parse_bool, default='False',
- help='Turn off BSM terms'
- )
- parser.add_argument(
- '--fix-mixing', type=misc_utils.parse_bool, default='False',
- help='Fix all mixing parameters to bi-maximal mixing'
- )
- parser.add_argument(
- '--fix-mixing-almost', type=misc_utils.parse_bool, default='False',
- help='Fix all mixing parameters except s23'
- )
- parser.add_argument(
- '--fix-scale', type=misc_utils.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(
- '--dimension', type=int, default=3,
- help='Set the new physics dimension to consider'
- )
- parser.add_argument(
- '--energy', type=float, default=1000,
- help='Set the energy scale'
- )
- parser.add_argument(
- '--seed', type=int, default=99,
- help='Set the random seed value'
- )
- parser.add_argument(
- '--threads', type=misc_utils.thread_type, default='1',
- help='Set the number of threads to use (int or "max")'
- )
- parser.add_argument(
- '--outfile', type=str, default='./untitled',
- help='Path to output chains'
- )
- llh_utils.likelihood_argparse(parser)
- mcmc_utils.mcmc_argparse(parser)
- gf_utils.gf_argparse(parser)
- plot_utils.plot_argparse(parser)
- nuisance_argparse(parser)
- 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)
-
- np.random.seed(args.seed)
-
- asimov_paramset, mcmc_paramset = get_paramsets(args)
- outfile = misc_utils.gen_outfile_name(args)
- print '== {0:<25} = {1}'.format('outfile', outfile)
-
- if args.run_mcmc:
- if args.likelihood is Likelihood.GOLEMFIT:
- fitter = gf_utils.setup_fitter(args, asimov_paramset)
- else:
- fitter = None
-
- ln_prob = partial(
- llh_utils.ln_prob, args=args, fitter=fitter,
- asimov_paramset=asimov_paramset, mcmc_paramset=mcmc_paramset
- )
-
- ndim = len(mcmc_paramset)
- if args.mcmc_seed_type == MCMCSeedType.UNIFORM:
- p0 = mcmc_utils.flat_seed(
- mcmc_paramset, nwalkers=args.nwalkers
- )
- elif args.mcmc_seed_type == MCMCSeedType.GAUSSIAN:
- p0 = mcmc_utils.gaussian_seed(
- mcmc_paramset, nwalkers=args.nwalkers
- )
-
- samples = mcmc_utils.mcmc(
- p0 = p0,
- ln_prob = ln_prob,
- ndim = ndim,
- nwalkers = args.nwalkers,
- burnin = args.burnin,
- nsteps = args.nsteps,
- threads = 1
- # TODO(shivesh): broken because you cannot pickle a GolemFitPy object
- # threads = misc_utils.thread_factors(args.threads)[0]
- )
- mcmc_utils.save_chains(samples, outfile)
-
- plot_utils.chainer_plot(
- infile = outfile+'.npy',
- outfile = outfile[:5]+outfile[5:].replace('data', 'plots'),
- outformat = ['pdf'],
- args = args,
- mcmc_paramset = mcmc_paramset
- )
-
- out = args.bayes_output+'/fr_fr_evidence' + misc_utils.gen_identifier(args)
- scan_scales = np.linspace(
- np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.bayes_bins
- )
- if args.run_bayes_factor:
- from scipy.optimize import minimize
- def fit_flags(flag_dict):
- flags = gf.FitParametersFlag()
- for key in flag_dict.iterkeys():
- flags.__setattr__(key, flag_dict[key])
- return flags
-
- args.likelihood = Likelihood.GF_FREQ
- default_flags = {
- # False means it's not fixed in minimization
- 'astroFlavorAngle1' : True,
- 'astroFlavorAngle2' : True,
- 'astroENorm' : False,
- 'astroMuNorm' : False,
- 'astroTauNorm' : False,
- 'convNorm' : False,
- 'promptNorm' : False,
- 'muonNorm' : False,
- 'astroNorm' : False,
- 'astroParticleBalance' : True,
- 'astroDeltaGamma' : False,
- 'cutoffEnergy' : True,
- 'CRDeltaGamma' : True,
- 'piKRatio' : True,
- 'NeutrinoAntineutrinoRatio' : True,
- 'darkNorm' : True,
- 'domEfficiency' : True,
- 'holeiceForward' : True,
- 'anisotropyScale' : True,
- 'astroNormSec' : True,
- 'astroDeltaGammaSec' : True
- }
-
- mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True)
- scale = mcmc_paramset.from_tag(ParamTag.SCALE)
-
- if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ:
- fitter = gf_utils.setup_fitter(args, asimov_paramset)
- fitter.SetFitParametersFlag(fit_flags(default_flags))
- else: fitter = None
-
- data = np.zeros((args.bayes_bins, 2))
- for s_idx, sc in enumerate(scan_scales):
- print '== SCALE = {0:.0E}'.format(np.power(10, sc))
- def fn(scen):
- theta = map(float, scen) + [0, sc]
- try:
- llh = llh_utils.triangle_llh(
- theta=theta, args=args, asimov_paramset=asimov_paramset,
- mcmc_paramset=mcmc_paramset_freq, fitter=fitter
- )
- print 'llh', llh
- except:
- print 'forbidden'
- return 9e99
- return -llh
- res = minimize(fun=fn, x0=[0.1, 0.1, 0.1], method='L-BFGS-B',
- bounds=[(0, 1), (0, 1), (0, 1)])
- llh = fn(res.x)
- data[s_idx][0] = sc
- data[s_idx][1] = llh
-
- misc_utils.make_dir(out)
- np.save(out+'.npy', np.array(data))
-
- # dirname = os.path.dirname(out)
- # plot_utils.bayes_factor_plot(
- # dirname=dirname, outfile=out, outformat=['png'], args=args
- # )
-
- out = args.angles_lim_output+'/fr_anfr_evidence' + misc_utils.gen_identifier(args)
- if args.run_angles_limit:
- def fit_flags(flag_dict):
- flags = gf.FitParametersFlag()
- for key in flag_dict.iterkeys():
- flags.__setattr__(key, flag_dict[key])
- return flags
-
- args.likelihood = Likelihood.GF_FREQ
- 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' : True,
- 'piKRatio' : True,
- 'NeutrinoAntineutrinoRatio' : True,
- 'darkNorm' : True,
- 'domEfficiency' : True,
- 'holeiceForward' : True,
- 'anisotropyScale' : True,
- 'astroNormSec' : True,
- 'astroDeltaGammaSec' : True
- }
-
- mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True)
- scenarios = [
- [np.sin(np.pi/2.)**2, 0, 0, 0],
- [0, np.cos(np.pi/2.)**4, 0, 0],
- [0, 0, np.sin(np.pi/2.)**2, 0],
- ]
-
- if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ:
- fitter = gf_utils.setup_fitter(args, asimov_paramset)
- fitter.SetFitParametersFlag(fit_flags(default_flags))
- else: fitter = None
-
- data = np.zeros((len(scenarios), args.bayes_bins, 2))
- mm_angles = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES)
- sc_angles = mcmc_paramset_freq.from_tag(ParamTag.SCALE)[0]
- for idx, scen in enumerate(scenarios):
- scales, llhs = [], []
- for yidx, an in enumerate(mm_angles):
- an.value = scen[yidx]
- for s_idx, sc in enumerate(scan_scales):
- print '== SCALE = {0:.0E}'.format(np.power(10, sc))
- theta = scen + [sc]
- print 'theta', theta
- llh = llh_utils.triangle_llh(
- theta=theta, args=args, asimov_paramset=asimov_paramset,
- mcmc_paramset=mcmc_paramset_freq, fitter=fitter
- )
- print 'llh', llh
- scales.append(sc)
- llhs.append(llh)
-
- for i, d in enumerate(llhs):
- data[idx][i][0] = scales[i]
- data[idx][i][1] = d
-
- misc_utils.make_dir(out)
- print 'saving to {0}.npy'.format(out)
- np.save(out+'.npy', np.array(data))
-
- dirname = os.path.dirname(out)
- plot_utils.plot_BSM_angles_limit(
- dirname=dirname, outfile=outfile, outformat=['png'],
- args=args, bayesian=True
- )
-
- out = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args)
- if args.run_angles_correlation:
- def fit_flags(flag_dict):
- flags = gf.FitParametersFlag()
- for key in flag_dict.iterkeys():
- flags.__setattr__(key, flag_dict[key])
- return flags
-
- args.likelihood = Likelihood.GF_FREQ
- 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' : True,
- 'piKRatio' : True,
- 'NeutrinoAntineutrinoRatio' : True,
- 'darkNorm' : True,
- 'domEfficiency' : True,
- 'holeiceForward' : True,
- 'anisotropyScale' : True,
- 'astroNormSec' : True,
- 'astroDeltaGammaSec' : True
- }
-
- mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True)
-
- scenarios = [
- [1, 0, 0, 0],
- [0, 1, 0, 0],
- [0, 0, 1, 0],
- ]
- mm_angles = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES)
- sc_angles = mcmc_paramset_freq.from_tag(ParamTag.SCALE)[0]
-
- if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ:
- fitter = gf_utils.setup_fitter(args, asimov_paramset)
- fitter.SetFitParametersFlag(fit_flags(default_flags))
- else: fitter = None
-
- scan_angles = np.linspace(0, 1, args.bayes_bins)
-
- data = np.zeros((len(scenarios), args.bayes_bins, args.bayes_bins, 3))
- for idx, scen in enumerate(scenarios):
- for an in mm_angles:
- an.value = 0
- keep = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES)[idx]
-
- for s_idx, sc in enumerate(scan_scales):
- scales, angles, llhs = [], [], []
- for a_idx, an in enumerate(scan_angles):
- print '== SCALE = {0:.0E}'.format(np.power(10, sc))
- print '== ANGLE = {0:.2f}'.format(an)
- sc_angles.value = sc
- s2_2 = an
- if s_idx == 0: keep.value = np.sin(np.arcsin(np.sqrt(s2_2))/2.)**2
- elif s_idx == 1: keep.value = np.cos(np.arcsin(np.sqrt(s2_2))/2.)**4
- elif s_idx == 2: keep.value = np.sin(np.arcsin(np.sqrt(s2_2))/2.)**2
- theta = mcmc_paramset_freq.values
- print 'mcmc_paramset_freq', mcmc_paramset_freq
- try:
- llh = llh_utils.triangle_llh(
- theta=theta, args=args, asimov_paramset=asimov_paramset,
- mcmc_paramset=mcmc_paramset_freq, fitter=fitter
- )
- except AssertionError:
- print 'forbidden by unitarity'
- data[idx][s_idx][a_idx][0] = np.nan
- data[idx][s_idx][a_idx][1] = np.nan
- data[idx][s_idx][a_idx][2] = np.nan
- continue
- print 'llh', llh
- data[idx][s_idx][a_idx][0] = sc
- data[idx][s_idx][a_idx][1] = an
- data[idx][s_idx][a_idx][2] = llh
-
- misc_utils.make_dir(out)
- print 'saving to {0}.npy'.format(out)
- np.save(out+'.npy', np.array(data))
-
- print "DONE!"
-
-
-main.__doc__ = __doc__
-
-
-if __name__ == '__main__':
- main()
-
diff --git a/sens.py b/sens.py
index 50fe0c8..f352149 100755
--- a/sens.py
+++ b/sens.py
@@ -84,6 +84,10 @@ def process_args(args):
'--sens-run and --fix-scale cannot be used together'
)
+ if args.sens_eval_bin is not None and args.plot_statistic:
+ print 'Cannot make plot when specific scale bin is chosen'
+ args.plot_statistic = False
+
args.measured_ratio = fr_utils.normalise_fr(args.measured_ratio)
if args.fix_source_ratio:
args.source_ratio = fr_utils.normalise_fr(args.source_ratio)
@@ -160,6 +164,7 @@ def parse_args(args=None):
if args is None: return parser.parse_args()
else: return parser.parse_args(args.split())
+
def main():
args = parse_args()
process_args(args)
@@ -192,16 +197,19 @@ def main():
np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.sens_bins
)
- if args.sens_eval_bin is None:
- eval_dim = args.sens_bins
- else: eval_dim = 1
+ corr_angles_categ = [SensitivityCateg.CORR_ANGLE, SensitivityCateg.CORR_ONE_ANGLE]
+ fixed_angle_categ = [SensitivityCateg.FIXED_ANGLE, SensitivityCateg.FIXED_ONE_ANGLE]
- if args.run_method is SensitivityCateg.CORR_ANGLE:
- scan_angles = np.linspace(0+1e-9, 1-1e-9, eval_dim)
+ if args.run_method in corr_angles_categ:
+ scan_angles = np.linspace(0+1e-9, 1-1e-9, args.sens_bins)
else: scan_angles = np.array([0])
print 'scan_scales', scan_scales
print 'scan_angles', scan_angles
+ if args.sens_eval_bin is None:
+ eval_dim = args.sens_bins
+ else: eval_dim = 1
+
out = args.outfile+'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
+ misc_utils.gen_identifier(args)
if args.sens_run:
@@ -218,18 +226,17 @@ def main():
if args.run_method is SensitivityCateg.FULL:
statistic_arr = np.full((eval_dim, 2), np.nan)
- elif args.run_method is SensitivityCateg.FIXED_ANGLE:
+ elif args.run_method in fixed_angle_categ:
statistic_arr = np.full((len(st_paramset_arr), eval_dim, 2), np.nan)
- elif args.run_method is SensitivityCateg.CORR_ANGLE:
+ elif args.run_method in corr_angles_categ:
statistic_arr = np.full(
(len(st_paramset_arr), eval_dim, eval_dim, 3), np.nan
)
for idx_scen, sens_paramset in enumerate(st_paramset_arr):
print '|||| SCENARIO = {0}'.format(idx_scen)
- if args.run_method in [SensitivityCateg.FIXED_ONE_ANGLE, SensitivityCateg.FIXED_ANGLE]:
- if SensitivityCateg.FIXED_ANGLE:
- for x in mmangles: x.value = 0.+1e-9
+ if args.run_method in fixed_angle_categ:
+ for x in mmangles: x.value = 0.+1e-9
if idx_scen == 0 or idx_scen == 2:
mmangles[idx_scen].value = np.sin(np.pi/4., dtype=DTYPE)**2
"""s_12^2 or s_23^2"""
@@ -238,22 +245,24 @@ def main():
"""c_13^4"""
for idx_an, an in enumerate(scan_angles):
- if args.run_method in [SensitivityCateg.CORR_ANGLE,
- SensitivityCateg.CORR_ONE_ANGLE]:
- print '|||| ANGLE = {0:<04.2}'.format(an)
- if SensitivityCateg.CORR_ANGLE:
- for x in mmangles: x.value = 0.+1e-9
+ if args.run_method in corr_angles_categ:
+ for x in mmangles: x.value = 0.+1e-9
mmangles[idx_scen].value = an
for idx_sc, sc in enumerate(scan_scales):
if args.sens_eval_bin is not None:
- if idx_sc == args.sens_eval_bin:
- out += '_scale_{0:.0E}'.format(np.power(10, sc))
- if args.run_method is SensitivityCateg.CORR_ANGLE:
- out += '_angle_{0:>03}'.format(int(an*100))
- break
+ if args.run_method in corr_angles_categ:
+ index = idx_an*args.sens_bins + idx_sc
+ else: index = idx_sc
+ if index == args.sens_eval_bin:
+ if idx_scen == 0:
+ out += '_scale_{0:.0E}'.format(np.power(10, sc))
+ if args.run_method in corr_angles_categ:
+ out += '_angle_{0:<04.2}'.format(an)
else: continue
+ if idx_sc == 0 or args.sens_eval_bin is not None:
+ print '|||| ANGLE = {0:<04.2}'.format(float(an))
print '|||| SCALE = {0:.0E}'.format(np.power(10, sc))
scale.value = sc
if args.stat_method is StatCateg.BAYESIAN:
@@ -296,10 +305,16 @@ def main():
print '=== final llh', stat
if args.run_method is SensitivityCateg.FULL:
statistic_arr[idx_sc] = np.array([sc, stat])
- elif args.run_method is SensitivityCateg.FIXED_ANGLE:
- statistic_arr[idx_scen][idx_sc] = np.array([sc, stat])
- elif args.run_method is SensitivityCateg.CORR_ANGLE:
- statistic_arr[idx_scen][idx_an][idx_sc] = np.array([an, sc, stat])
+ elif args.run_method in fixed_angle_categ:
+ if args.sens_eval_bin is not None:
+ statistic_arr[idx_scen][0] = np.array([sc, stat])
+ else:
+ statistic_arr[idx_scen][idx_sc] = np.array([sc, stat])
+ elif args.run_method in corr_angles_categ:
+ if args.sens_eval_bin is not None:
+ statistic_arr[idx_scen][0][0] = np.array([an, sc, stat])
+ else:
+ statistic_arr[idx_scen][idx_an][idx_sc] = np.array([an, sc, stat])
misc_utils.make_dir(out)
print 'Saving to {0}'.format(out+'.npy')
@@ -310,7 +325,7 @@ def main():
else: raw = np.load(out+'.npy')
data = ma.masked_invalid(raw, 0)
- basename = os.path.dirname(out) + '/mnrun/' + os.path.basename(out)
+ basename = os.path.dirname(out) + '/statrun/' + os.path.basename(out)
baseoutfile = basename[:5]+basename[5:].replace('data', 'plots')
if args.run_method is SensitivityCateg.FULL:
plot_utils.plot_statistic(
@@ -320,7 +335,7 @@ def main():
args = args,
scale_param = scale
)
- elif args.run_method is SensitivityCateg.FIXED_ANGLE:
+ elif args.run_method in fixed_angle_categ:
for idx_scen in xrange(len(st_paramset_arr)):
print '|||| SCENARIO = {0}'.format(idx_scen)
outfile = baseoutfile + '_SCEN{0}'.format(idx_scen)
@@ -335,7 +350,7 @@ def main():
scale_param = scale,
label = label
)
- elif args.run_method is SensitivityCateg.CORR_ANGLE:
+ elif args.run_method in corr_angles_categ:
for idx_scen in xrange(len(st_paramset_arr)):
print '|||| SCENARIO = {0}'.format(idx_scen)
basescenoutfile = baseoutfile + '_SCEN{0}'.format(idx_scen)
@@ -343,11 +358,11 @@ def main():
elif idx_scen == 1: label = r'$\mathcal{O}_{13}='
elif idx_scen == 2: label = r'$\mathcal{O}_{23}='
for idx_an, an in enumerate(scan_angles):
- print '|||| ANGLE = {0:<04.2}'.format(an)
+ print '|||| ANGLE = {0:<04.2}'.format(float(an))
outfile = basescenoutfile + '_ANGLE{0}'.format(idx_an)
_label = label + r'{0:<04.2}$'.format(an)
plot_utils.plot_statistic(
- data = data[idx_scen][idx_an][:,1:],
+ data = data[idx_scen][index][:,1:],
outfile = outfile,
outformat = ['png'],
args = args,
@@ -361,4 +376,3 @@ main.__doc__ = __doc__
if __name__ == '__main__':
main()
-
diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py
index 79f9b6d..2866887 100644
--- a/submitter/mcmc_dag.py
+++ b/submitter/mcmc_dag.py
@@ -9,9 +9,9 @@ full_scan_mfr = [
fix_sfr_mfr = [
(1, 1, 1, 1, 2, 0),
- (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, 0, 1, 2, 0),
# (1, 1, 0, 1, 0, 0),
# (1, 1, 0, 0, 1, 0),
@@ -26,10 +26,10 @@ GLOBAL_PARAMS = {}
# MCMC
GLOBAL_PARAMS.update(dict(
run_mcmc = 'True',
- burnin = 2500,
- nsteps = 10000,
+ burnin = 250,
+ nsteps = 1000,
nwalkers = 60,
- seed = 'None',
+ seed = 25,
mcmc_seed_type = 'uniform'
))
@@ -90,9 +90,9 @@ with open(outfile, 'w') as f:
f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, frs[3]))
f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, frs[4]))
f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, frs[5]))
- for key in GLOBAL_PARAMS.keys():
- f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, GLOBAL_PARAMS[key]))
- f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outfile))
+ for key in GLOBAL_PARAMS.iterkeys():
+ f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key]))
+ f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
job_number += 1
for frs in full_scan_mfr:
@@ -110,7 +110,7 @@ with open(outfile, 'w') as f:
f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, 0))
f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, 0))
f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, 0))
- for key in GLOBAL_PARAMS.keys():
- f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, GLOBAL_PARAMS[key]))
- f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outfile))
+ for key in GLOBAL_PARAMS.iterkeys():
+ f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key]))
+ f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
job_number += 1
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index 2cecfe9..63aaaba 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -9,9 +9,9 @@ full_scan_mfr = [
fix_sfr_mfr = [
(1, 1, 1, 1, 2, 0),
- (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, 0, 1, 2, 0),
# (1, 1, 0, 1, 0, 0),
# (1, 1, 0, 0, 1, 0),
@@ -24,12 +24,13 @@ fix_sfr_mfr = [
GLOBAL_PARAMS = {}
# Bayes Factor
+sens_eval_bin = 'all' # set to 'all' to run normally
GLOBAL_PARAMS.update(dict(
sens_run = 'True',
- run_method = 'full',
- stat_method = 'bayesian',
+ run_method = 'corr_angle',
+ stat_method = 'frequentist',
sens_bins = 10,
- sens_eval_bin = 'all' # set to 'all' to run normally
+ seed = 'None'
))
# MultiNest
@@ -69,11 +70,13 @@ GLOBAL_PARAMS.update(dict(
plot_statistic = 'True'
))
-outfile = 'dagman_FR_SENS.submit'
+outfile = 'dagman_FR_SENS_{0}_{1}.submit'.format(
+ GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method']
+)
golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
-if GLOBAL_PARAMS['sens_eval_bin'].lower() != 'all':
+if sens_eval_bin.lower() != 'all':
if GLOBAL_PARAMS['run_method'].lower() == 'corr_angle':
sens_runs = GLOBAL_PARAMS['sens_bins']**2
else:
@@ -105,9 +108,9 @@ with open(outfile, 'w') as f:
f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, frs[4]))
f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, frs[5]))
f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r))
- for key in GLOBAL_PARAMS.keys():
- f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, GLOBAL_PARAMS[key]))
- f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outfile))
+ for key in GLOBAL_PARAMS.iterkeys():
+ f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key]))
+ f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, output))
job_number += 1
for frs in full_scan_mfr:
@@ -128,7 +131,7 @@ with open(outfile, 'w') as f:
f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, 0))
f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, 0))
f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r))
- for key in GLOBAL_PARAMS.keys():
- f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, GLOBAL_PARAMS[key]))
- f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outfile))
+ for key in GLOBAL_PARAMS.iterkeys():
+ f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key]))
+ f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, output))
job_number += 1
diff --git a/submitter/sens_submit.sub b/submitter/sens_submit.sub
index 1a02608..57a2588 100644
--- a/submitter/sens_submit.sub
+++ b/submitter/sens_submit.sub
@@ -1,4 +1,4 @@
-Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/fr.py
+Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/sens.py
Arguments = "--ast $(ast) --data $(data) --dimension $(dimension) --fix-mixing $(fix_mixing) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --no-bsm $(no_bsm) --outfile $(outfile) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning) --fix-mixing-almost $(fix_mixing_almost) --sens-run $(sens_run) --run-method $(run_method) --stat-method $(stat_method) --sens-bins $(sens_bins) --sens-eval-bin $(sens_eval_bin) --mn-live-points $(mn_live_points) --mn-tolerance $(mn_tolerance) --mn-output $(mn_output) --plot-statistic $(plot_statistic)"
# All logs will go to a single file
diff --git a/utils/fr.py b/utils/fr.py
index a82e081..b2a1274 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -9,7 +9,6 @@ Useful functions for the BSM flavour ratio analysis
from __future__ import absolute_import, division
-import argparse
from functools import partial
import numpy as np
@@ -354,7 +353,7 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
if np.shape(sm_u) != (3, 3):
raise ValueError(
'Input matrix should be a square and dimension 3, '
- 'got\n{0}'.format(ham)
+ 'got\n{0}'.format(sm_u)
)
if fix_mixing and fix_mixing_almost:
diff --git a/utils/gf.py b/utils/gf.py
index 17ac029..0770401 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -9,9 +9,10 @@ Useful GolemFit wrappers for the BSM flavour ratio analysis
from __future__ import absolute_import, division
-import socket
from functools import partial
+import numpy as np
+
import GolemFitPy as gf
from utils.enums import DataType, SteeringCateg
diff --git a/utils/likelihood.py b/utils/likelihood.py
index 6387a1e..9629b65 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -14,8 +14,6 @@ from functools import partial
import numpy as np
from scipy.stats import multivariate_normal, rv_continuous
-import GolemFitPy as gf
-
from utils import fr as fr_utils
from utils import gf as gf_utils
from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg
@@ -89,6 +87,10 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
if args.energy_dependance is EnergyDependance.SPECTRAL:
bin_centers = np.sqrt(args.binning[:-1]*args.binning[1:])
bin_width = np.abs(np.diff(args.binning))
+ if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
+ if 'astroDeltaGamma' in hypo_paramset.names:
+ args.spectral_index = hypo_paramset['astroDeltaGamma'].value
+ print 'args.spectral_index', args.spectral_index
if args.fix_source_ratio:
if args.energy_dependance is EnergyDependance.MONO:
diff --git a/utils/misc.py b/utils/misc.py
index 970c693..cad03bc 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -14,7 +14,6 @@ import errno
import multiprocessing
import argparse
-from collections import Sequence
from operator import attrgetter
import numpy as np
diff --git a/utils/multinest.py b/utils/multinest.py
index 005a43a..9dd0742 100644
--- a/utils/multinest.py
+++ b/utils/multinest.py
@@ -16,7 +16,7 @@ import numpy as np
from pymultinest import analyse, run
from utils import likelihood
-from utils.misc import gen_outfile_name, make_dir, parse_bool
+from utils.misc import gen_outfile_name, make_dir
def CubePrior(cube, ndim, n_params):
diff --git a/utils/plot.py b/utils/plot.py
index 0c431cf..0160da4 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -10,7 +10,6 @@ Plotting functions for the BSM flavour ratio analysis
from __future__ import absolute_import, division
import os
-import argparse
import numpy as np
import matplotlib as mpl
@@ -19,9 +18,7 @@ from matplotlib import rc
from matplotlib import pyplot as plt
from matplotlib.offsetbox import AnchoredText
-import getdist
-from getdist import plots
-from getdist import mcsamples
+from getdist import plots, mcsamples
from utils import misc as misc_utils
from utils.enums import EnergyDependance, Likelihood, ParamTag, StatCateg
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-17 09:50:31 +0000