From c37932036698600c7b44d2ff15aac6784d201098 Mon Sep 17 00:00:00 2001 From: shivesh Date: Sat, 28 Apr 2018 17:01:52 -0500 Subject: Sat Apr 28 17:01:52 CDT 2018 --- fr.py | 5 +- fr_edit.py | 581 ---------------------------------------------- sens.py | 74 +++--- submitter/mcmc_dag.py | 24 +- submitter/sens_dag.py | 31 +-- submitter/sens_submit.sub | 2 +- utils/fr.py | 3 +- utils/gf.py | 3 +- utils/likelihood.py | 6 +- utils/misc.py | 1 - utils/multinest.py | 2 +- utils/plot.py | 5 +- 12 files changed, 87 insertions(+), 650 deletions(-) delete mode 100755 fr_edit.py 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 -- cgit v1.2.3