diff options
| -rwxr-xr-x | fr.py | 42 | ||||
| -rwxr-xr-x | plot_sens.py | 368 | ||||
| -rwxr-xr-x | sens.py | 46 | ||||
| -rw-r--r-- | submitter/mcmc_dag.py | 7 | ||||
| -rw-r--r-- | submitter/mcmc_submit.sub | 6 | ||||
| -rw-r--r-- | submitter/sens_dag.py | 21 | ||||
| -rw-r--r-- | submitter/sens_submit.sub | 4 | ||||
| -rw-r--r-- | test/test_LV.py | 96 | ||||
| -rw-r--r-- | test/test_NSI.py | 106 | ||||
| -rw-r--r-- | test/test_all_gf.py | 79 | ||||
| -rw-r--r-- | test/test_gf.py | 21 | ||||
| -rw-r--r-- | test/test_gf_simple.py | 40 | ||||
| -rw-r--r-- | utils/fr.py | 2 | ||||
| -rw-r--r-- | utils/gf.py | 5 | ||||
| -rw-r--r-- | utils/likelihood.py | 4 | ||||
| -rw-r--r-- | utils/misc.py | 30 | ||||
| -rw-r--r-- | utils/multinest.py | 10 | ||||
| -rw-r--r-- | utils/param.py | 14 | ||||
| -rw-r--r-- | utils/plot.py | 381 |
19 files changed, 1111 insertions, 171 deletions
@@ -33,28 +33,28 @@ def define_nuisance(): g_prior = PriorsCateg.GAUSSIAN hg_prior = PriorsCateg.HALFGAUSS e = 1e-9 - # nuisance.extend([ - # Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), - # Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.995, 1-e], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=hg_prior, tag=tag), - # Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), - # Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), - # Param( - # name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23], - # std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag - # ), - # Param( - # name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21], - # std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag - # ) - # ]) + nuisance.extend([ + Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), + Param(name='c_13_4', value=(1-(0.02206))**2, seed=[0.950, 0.961], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=g_prior, tag=tag), + Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), + Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), + Param( + name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23], + std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag + ), + Param( + name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21], + std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag + ) + ]) tag = ParamTag.NUISANCE - # nuisance.extend([ - # Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3, tag=tag), - # Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag), - # Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1, tag=tag), - # Param(name='astroNorm', value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1, tag=tag), - # Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) - # ]) + nuisance.extend([ + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + ]) return ParamSet(nuisance) diff --git a/plot_sens.py b/plot_sens.py index e51d55d..0c862eb 100755 --- a/plot_sens.py +++ b/plot_sens.py @@ -10,21 +10,355 @@ HESE BSM flavour ratio analysis plotting script from __future__ import absolute_import, division +import os +import argparse +from functools import partial +from copy import deepcopy + import numpy as np -import matplotlib as mpl -mpl.use('Agg') -from matplotlib import rc -from matplotlib import pyplot as plt -from matplotlib.offsetbox import AnchoredText - -rc('text', usetex=False) -rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) - -FRS = [ - (1, 1, 1, 1, 2, 0), - (1, 1, 1, 0, 1, 0), - # (1, 1, 1, 1, 0, 0), - # (1, 1, 1, 0, 0, 1), -] - -DIMENSION = [3, 6] +import numpy.ma as ma + +from utils import fr as fr_utils +from utils import likelihood as llh_utils +from utils import misc as misc_utils +from utils import plot as plot_utils +from utils.enums import EnergyDependance, Likelihood, ParamTag +from utils.enums import PriorsCateg, SensitivityCateg, StatCateg +from utils.param import Param, ParamSet, get_paramsets + +from utils import multinest as mn_utils + + +def define_nuisance(): + """Define the nuisance parameters.""" + tag = ParamTag.SM_ANGLES + g_prior = PriorsCateg.GAUSSIAN + hg_prior = PriorsCateg.HALFGAUSS + e = 1e-9 + nuisance = [ + Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), + Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.995, 1-e], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=hg_prior, tag=tag), + Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), + Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), + Param( + name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23], + std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag + ), + Param( + name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21], + std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag + ) + ] + tag = ParamTag.NUISANCE + nuisance.extend([ + Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3, tag=tag), + Param(name='promptNorm', value=0., seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag), + Param(name='muonNorm', value=1., seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1, tag=tag), + Param(name='astroNorm', value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1, tag=tag), + Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1, tag=tag) + ]) + return ParamSet(nuisance) + + +def nuisance_argparse(parser): + nuisance = define_nuisance() + for parm in nuisance: + parser.add_argument( + '--'+parm.name, type=float, default=parm.value, + help=parm.name+' to inject' + ) + + +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.fix_scale: + raise NotImplementedError( + '--fix-scale not implemented' + ) + + args.measured_ratio = fr_utils.normalise_fr(args.measured_ratio) + if args.fix_source_ratio: + assert len(args.source_ratios) % 3 == 0 + srs = [args.source_ratios[3*x:3*x+3] + for x in range(int(len(args.source_ratios)/3))] + args.source_ratios = map(fr_utils.normalise_fr, srs) + + 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 args.split_jobs and args.run_method is SensitivityCateg.FULL: + raise NotImplementedError( + 'split_jobs and run_method not implemented' + ) + + args.dimensions = np.sort(args.dimensions) + + args_copy = deepcopy(args) + scale_regions = [] + for dim in args.dimensions: + args_copy.dimension = dim + _, scale_region = fr_utils.estimate_scale(args_copy) + scale_regions.append(scale_region) + args.scale_region = [np.min(scale_regions), np.max(scale_regions)] + args.scale = np.power(10., np.average(np.log10(args.scale_region))) + + +def parse_args(args=None): + """Parse command line arguments""" + parser = argparse.ArgumentParser( + description="HESE BSM flavour ratio analysis plotting script", + formatter_class=misc_utils.SortingHelpFormatter, + ) + parser.add_argument( + '--infile', type=str, default='./untitled', + help='Path to input dir' + ) + parser.add_argument( + '--run-method', default='full', + type=partial(misc_utils.enum_parse, c=SensitivityCateg), + choices=SensitivityCateg, + help='Choose which type of sensivity plot to make' + ) + parser.add_argument( + '--stat-method', default='bayesian', + type=partial(misc_utils.enum_parse, c=StatCateg), choices=StatCateg, + help='Statistical method to employ' + ) + parser.add_argument( + '--sens-bins', type=int, default=10, + help='Number of bins for the Bayes factor plot' + ) + parser.add_argument( + '--split-jobs', type=misc_utils.parse_bool, default='False', + help='Did the jobs get split' + ) + parser.add_argument( + '--plot', type=misc_utils.parse_bool, default='True', + help='Make sensitivity plots' + ) + parser.add_argument( + '--plot-statistic', type=misc_utils.parse_bool, default='False', + help='Plot MultiNest evidence or LLH value' + ) + fr_utils.fr_argparse(parser) + llh_utils.likelihood_argparse(parser) + mn_utils.mn_argparse(parser) + nuisance_argparse(parser) + misc_utils.remove_option(parser, 'dimension') + misc_utils.remove_option(parser, 'source_ratio') + misc_utils.remove_option(parser, 'scale') + misc_utils.remove_option(parser, 'scale_region') + parser.add_argument( + '--dimensions', type=int, nargs='*', default=[3, 6], + help='Set the new physics dimensions to consider' + ) + parser.add_argument( + '--source-ratios', type=int, nargs='*', default=[2, 1, 0], + help='Set the source flavour ratios for the case when you want to fix it' + ) + 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) + + asimov_paramset, llh_paramset = get_paramsets(args, define_nuisance()) + + scale = llh_paramset.from_tag(ParamTag.SCALE)[0] + mmangles = llh_paramset.from_tag(ParamTag.MMANGLES) + if args.run_method is SensitivityCateg.FULL: + st_paramset_arr = [llh_paramset.from_tag(ParamTag.SCALE, invert=True)] + elif args.run_method in [SensitivityCateg.FIXED_ANGLE, SensitivityCateg.CORR_ANGLE]: + nscale_pmset = llh_paramset.from_tag([ParamTag.SCALE, ParamTag.MMANGLES], invert=True) + st_paramset_arr = [nscale_pmset] * 3 + elif args.run_method in [SensitivityCateg.FIXED_ONE_ANGLE, SensitivityCateg.CORR_ONE_ANGLE]: + nscale_pmset = llh_paramset.from_tag(ParamTag.SCALE, invert=True) + st_paramset_arr = [] + for x in xrange(3): + st_paramset_arr.append( + ParamSet([prms for prms in nscale_pmset + if mmangles[x].name != prms.name]) + ) + + corr_angles_categ = [SensitivityCateg.CORR_ANGLE, SensitivityCateg.CORR_ONE_ANGLE] + fixed_angle_categ = [SensitivityCateg.FIXED_ANGLE, SensitivityCateg.FIXED_ONE_ANGLE] + + 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_angles', scan_angles + + dims = len(args.dimensions) + srcs = len(args.source_ratios) + if args.run_method is SensitivityCateg.FULL: + statistic_arr = np.full((dims, srcs, args.sens_bins, 2), np.nan) + elif args.run_method in fixed_angle_categ: + statistic_arr = np.full((dims, srcs, len(st_paramset_arr), args.sens_bins, 2), np.nan) + elif args.run_method in corr_angles_categ: + statistic_arr = np.full( + (dims, srcs, len(st_paramset_arr), args.sens_bins, args.sens_bins, 3), np.nan + ) + + print 'Loading data' + for idim, dim in enumerate(args.dimensions): + argsc = deepcopy(args) + argsc.dimension = dim + _, scale_region = fr_utils.estimate_scale(argsc) + argsc.scale_region = scale_region + scan_scales = np.linspace( + np.log10(scale_region[0]), np.log10(scale_region[1]), args.sens_bins + ) + + for isrc, src in enumerate(args.source_ratios): + argsc.source_ratio = src + infile = args.infile + if args.likelihood is Likelihood.GOLEMFIT: + infile += '/golemfit/' + elif args.likelihood is Likelihood.GAUSSIAN: + infile += '/gaussian/' + infile += '/DIM{0}/fix_ifr/'.format(dim) + if args.likelihood is Likelihood.GAUSSIAN: + infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_')) + infile += 'fr_stat/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \ + + misc_utils.gen_identifier(argsc) + print '== {0:<25} = {1}'.format('infile', infile) + + if args.split_jobs: + for idx_an, an in enumerate(scan_angles): + for idx_sc, sc in enumerate(scan_scales): + filename = infile + '_scale_{0:.0E}'.format(np.power(10, sc)) + try: + if args.run_method in fixed_angle_categ: + print 'Loading from {0}'.format(filename+'.npy') + statistic_arr[idim][isrc][:,idx_sc] = np.load(filename+'.npy')[:,0] + if args.run_method in corr_angles_categ: + filename += '_angle_{0:<04.2}'.format(an) + print 'Loading from {0}'.format(filename+'.npy') + statistic_arr[idim][isrc][:,idx_an,idx_sc] = np.load(filename+'.npy')[:,0,0] + except: + print 'Unable to load file {0}'.format(filename+'.npy') + continue + else: + print 'Loading from {0}'.format(infile+'.npy') + try: + statistic_arr[idim][isrc] = np.load(infile+'.npy') + except: + print 'Unable to load file {0}'.format(infile+'.npy') + continue + + print 'statistic_arr', statistic_arr + + data = ma.masked_invalid(statistic_arr) + if args.plot_statistic: + print 'Plotting statistic' + + argsc = deepcopy(args) + base_infile = args.infile + if args.likelihood is Likelihood.GOLEMFIT: + base_infile += '/golemfit/' + elif args.likelihood is Likelihood.GAUSSIAN: + base_infile += '/gaussian/' + for idim, dim in enumerate(args.dimensions): + argsc.dimension = dim + _, scale_region = fr_utils.estimate_scale(argsc) + argsc.scale_region = scale_region + infile = base_infile + '/DIM{0}/fix_ifr/'.format(dim) + if args.likelihood is Likelihood.GAUSSIAN: + infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_')) + infile += 'fr_stat/' + + for isrc, src in enumerate(args.source_ratios): + argsc.source_ratio = src + finfile = infile +'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \ + + misc_utils.gen_identifier(argsc) + basename = os.path.dirname(finfile) + '/statrun/' + os.path.basename(finfile) + baseoutfile = basename[:5]+basename[5:].replace('data', 'plots') + if args.run_method is SensitivityCateg.FULL: + outfile = baseoutfile + plot_utils.plot_statistic( + data = data[idim][isrc], + outfile = outfile, + outformat = ['png'], + args = argsc, + scale_param = scale, + ) + if 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) + if idx_scen == 0: label = r'$\mathcal{O}_{12}=\frac{\pi}{4}$' + elif idx_scen == 1: label = r'$\mathcal{O}_{13}=\frac{\pi}{4}$' + elif idx_scen == 2: label = r'$\mathcal{O}_{23}=\frac{\pi}{4}$' + plot_utils.plot_statistic( + data = data[idim][isrc][idx_scen], + outfile = outfile, + outformat = ['png'], + args = argsc, + scale_param = scale, + label = label + ) + 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) + if idx_scen == 0: label = r'$\mathcal{O}_{12}=' + 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(float(an)) + outfile = basescenoutfile + '_ANGLE{0}'.format(idx_an) + _label = label + r'{0:<04.2}$'.format(an) + plot_utils.plot_statistic( + data = data[idim][isrc][idx_scen][idx_an][:,1:], + outfile = outfile, + outformat = ['png'], + args = argsc, + scale_param = scale, + label = _label + ) + + if args.plot: + print 'Plotting sensitivities' + + basename = args.infile + '{0}/{1}'.format(args.likelihood, args.stat_method) + baseoutfile = basename[:5]+basename[5:].replace('data', 'plots') + + if args.run_method is SensitivityCateg.FULL: + plot_utils.plot_sens_full( + data = data, + outfile = basename + '/FULL', + outformat = ['png'], + args = args, + ) + elif args.run_method in fixed_angle_categ: + plot_utils.plot_sens_fixed_angle( + data = data, + outfile = basename + '/FIXED_ANGLE', + outformat = ['png'], + args = args, + ) + elif args.run_method in corr_angles_categ: + plot_utils.plot_sens_corr_angle( + data = data, + outfile = basename + '/CORR_ANGLE', + outformat = ['png'], + args = args, + ) + + +main.__doc__ = __doc__ + + +if __name__ == '__main__': + main() @@ -34,14 +34,14 @@ from utils import multinest as mn_utils def define_nuisance(): """Define the nuisance parameters.""" tag = ParamTag.SM_ANGLES + nuisance = [] g_prior = PriorsCateg.GAUSSIAN - hg_prior = PriorsCateg.HALFGAUSS e = 1e-9 - nuisance = [ - Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), - Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.995, 1-e], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=hg_prior, tag=tag), - Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), - Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), + nuisance.extend([ + Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), + Param(name='c_13_4', value=(1-(0.02206))**2, seed=[0.950, 0.961], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=g_prior, tag=tag), + Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), + Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), Param( name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23], std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag @@ -50,7 +50,7 @@ def define_nuisance(): name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21], std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag ) - ] + ]) tag = ParamTag.NUISANCE nuisance.extend([ Param(name='convNorm', value=1., seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3, tag=tag), @@ -173,10 +173,8 @@ def main(): np.random.seed(args.seed) asimov_paramset, llh_paramset = get_paramsets(args, define_nuisance()) - scale = llh_paramset.from_tag(ParamTag.SCALE)[0] - outfile = misc_utils.gen_outfile_name(args) - print '== {0:<25} = {1}'.format('outfile', outfile) + scale = llh_paramset.from_tag(ParamTag.SCALE)[0] mmangles = llh_paramset.from_tag(ParamTag.MMANGLES) if args.run_method is SensitivityCateg.FULL: st_paramset_arr = [llh_paramset.from_tag(ParamTag.SCALE, invert=True)] @@ -246,7 +244,13 @@ def main(): for idx_an, an in enumerate(scan_angles): if args.run_method in corr_angles_categ: for x in mmangles: x.value = 0.+1e-9 - mmangles[idx_scen].value = an + angle = np.arcsin(np.sqrt(an))/2. + if idx_scen == 0 or idx_scen == 2: + mmangles[idx_scen].value = np.sin(angle)**2 + """s_12^2 or s_23^2""" + elif idx_scen == 1: + mmangles[idx_scen].value = np.cos(angle)**4 + """c_13^4""" for idx_sc, sc in enumerate(scan_scales): if args.sens_eval_bin is not None: @@ -275,6 +279,7 @@ def main(): ) except: print 'Failed run, continuing' + # raise continue print '## Evidence = {0}'.format(stat) elif args.stat_method is StatCateg.FREQUENTIST: @@ -289,18 +294,27 @@ def main(): llh_paramset=llh_paramset, fitter=fitter ) # print 'llh_paramset', llh_paramset - # print 'llh', llh + # print '-llh', -llh return -llh n_params = len(sens_paramset) - x0 = np.full(n_params, 0.7) bounds = np.dstack([np.zeros(n_params), np.ones(n_params)])[0] + x0_arr = np.linspace(0+1e-3, 1-1e-3, 3) + stat = -np.inf try: - res = minimize(fun=fn, x0=x0, method='L-BFGS-B', bounds=bounds) + for x0_v in x0_arr: + x0 = np.full(n_params, x0_v) + # res = minimize(fun=fn, x0=x0, method='Powell', tol=1e-3) + res = minimize(fun=fn, x0=x0, method='Nelder-Mead', tol=1e-3) + # res = minimize(fun=fn, x0=x0, method='L-BFGS-B', tol=1e-3, bounds=bounds) + # res = minimize(fun=fn, x0=x0, method='BFGS', tol=1e-3) + s = -fn(res.x) + if s > stat: + stat = s except AssertionError: print 'Failed run, continuing' + # raise continue - stat = -fn(res.x) print '=== final llh', stat if args.run_method is SensitivityCateg.FULL: statistic_arr[idx_sc] = np.array([sc, stat]) @@ -323,7 +337,7 @@ def main(): print 'Plotting statistic' if args.sens_run: raw = statistic_arr else: raw = np.load(out+'.npy') - data = ma.masked_invalid(raw, 0) + data = ma.masked_invalid(raw) basename = os.path.dirname(out) + '/statrun/' + os.path.basename(out) baseoutfile = basename[:5]+basename[5:].replace('data', 'plots') diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py index 7f6eae6..d06e920 100644 --- a/submitter/mcmc_dag.py +++ b/submitter/mcmc_dag.py @@ -35,6 +35,7 @@ GLOBAL_PARAMS.update(dict( # FR dimension = [3, 6] +# dimension = [4, 5, 7, 8] GLOBAL_PARAMS.update(dict( threads = 1, binning = '1e4 1e7 5', @@ -44,7 +45,7 @@ GLOBAL_PARAMS.update(dict( spectral_index = -2, fix_mixing = 'False', fix_mixing_almost = 'False', - fold_index = 'False' + fold_index = 'True' )) # Likelihood @@ -65,7 +66,7 @@ GLOBAL_PARAMS.update(dict( plot_elements = 'False', )) -outfile = 'dagman_FR_MCMC.submit' +outfile = 'dagman_FR_MCMC_{0}.submit'.format(GLOBAL_PARAMS['likelihood']) golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit' condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/mcmc_submit.sub' @@ -115,3 +116,5 @@ with open(outfile, 'w') as f: 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 + + print 'dag file = {0}'.format(outfile) diff --git a/submitter/mcmc_submit.sub b/submitter/mcmc_submit.sub index f006350..c565205 100644 --- a/submitter/mcmc_submit.sub +++ b/submitter/mcmc_submit.sub @@ -1,5 +1,5 @@ Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/fr.py -Arguments = "--ast $(ast) --burnin $(burnin) --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) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --run-mcmc $(run_mcmc) --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) --fold-index $(fold_index)" +Arguments = "--ast $(ast) --burnin $(burnin) --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) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --run-mcmc $(run_mcmc) --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) --fold-index $(fold_index) --mcmc-seed-type $(mcmc_seed_type)" # All logs will go to a single file log = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log @@ -16,7 +16,7 @@ getenv = True # +TransferOutput="" Transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/ -request_memory = 2GB +request_memory = 3GB request_cpus = 1 Universe = vanilla @@ -33,7 +33,7 @@ Notification = never # Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6") # Requirements = IS_GLIDEIN -# Requirements = (OpSysMajorVer =?= 6) +Requirements = (OpSysMajorVer =?= 6) # GO! queue diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py index 17b063d..00be770 100644 --- a/submitter/sens_dag.py +++ b/submitter/sens_dag.py @@ -27,21 +27,23 @@ GLOBAL_PARAMS = {} sens_eval_bin = 'all' # set to 'all' to run normally GLOBAL_PARAMS.update(dict( sens_run = 'True', - run_method = 'fixed_angle', # full, fixed_angle, corr_angle + run_method = 'full', # full, fixed_angle, corr_angle stat_method = 'frequentist', - sens_bins = 10, + sens_bins = 40, seed = 'None' )) # MultiNest GLOBAL_PARAMS.update(dict( - mn_live_points = 400, + mn_live_points = 800, mn_tolerance = 0.01, mn_output = './mnrun' )) # FR -dimension = [3, 6] +# dimension = [3, 6] +dimension = [4, 5, 7, 8] +# dimension = [3, 4, 5, 6, 7, 8] GLOBAL_PARAMS.update(dict( threads = 1, binning = '1e4 1e7 5', @@ -50,12 +52,13 @@ GLOBAL_PARAMS.update(dict( energy_dependance = 'spectral', spectral_index = -2, fix_mixing = 'False', - fix_mixing_almost = 'False' + fix_mixing_almost = 'False', + fold_index = 'False' )) # Likelihood GLOBAL_PARAMS.update(dict( - likelihood = 'golemfit', + likelihood = 'gaussian', sigma_ratio = '0.01' )) @@ -70,8 +73,8 @@ GLOBAL_PARAMS.update(dict( plot_statistic = 'True' )) -outfile = 'dagman_FR_SENS_{0}_{1}.submit'.format( - GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'] +outfile = 'dagman_FR_SENS_{0}_{1}_{2}.submit'.format( + GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'], GLOBAL_PARAMS['likelihood'] ) golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit' condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub' @@ -141,3 +144,5 @@ with open(outfile, 'w') as f: 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 + + print 'dag file = {0}'.format(outfile) diff --git a/submitter/sens_submit.sub b/submitter/sens_submit.sub index 44ec0a3..f3a3e4c 100644 --- a/submitter/sens_submit.sub +++ b/submitter/sens_submit.sub @@ -16,7 +16,7 @@ getenv = True # +TransferOutput="" Transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/ -request_memory = 1GB +request_memory = 3GB request_cpus = 1 Universe = vanilla @@ -33,7 +33,7 @@ Notification = never # Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6") # Requirements = IS_GLIDEIN -# Requirements = (OpSysMajorVer =?= 6) +Requirements = (OpSysMajorVer =?= 6) # GO! queue diff --git a/test/test_LV.py b/test/test_LV.py new file mode 100644 index 0000000..72d0a9c --- /dev/null +++ b/test/test_LV.py @@ -0,0 +1,96 @@ +#!/usr/bin/env python + +from __future__ import absolute_import, division + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rc + +import GolemFitPy as gf + +rc('text', usetex=True) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +dp = gf.DataPaths() +steer = gf.SteeringParams() + +fig = plt.figure(figsize=[12, 8]) +ax = fig.add_subplot(111) +ax.set_xscale('log') +ax.set_yscale('log') + +npp = gf.NewPhysicsParams() +npp.type = gf.NewPhysicsType.None +# npp.n_lv = 1 +# npp.lambda_1 = 1.e100 +# npp.lambda_2 = 1.e100 + +golem = gf.GolemFit(dp, steer, npp) + +binning = golem.GetEnergyBinsMC() +ax.set_xlim(binning[0], binning[-1]) +# ax.set_ylim(binning[0], binning[-1]) + +fit_params = gf.FitParameters(gf.sampleTag.HESE) +golem.SetupAsimov(fit_params) + +exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, + drawstyle='steps-pre', label='NULL', linestyle='--') + +print 'NULL min_llh', golem.MinLLH().likelihood +print 'NULL expectation', exp +print + +npp = gf.NewPhysicsParams() +npp.type = gf.NewPhysicsType.LorentzViolation +npp.n_lv = 1 +npp.lambda_1 = 1.e20 +npp.lambda_2 = 1.e20 + +golem.SetNewPhysicsParams(npp) + +exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, + drawstyle='steps-pre', label='1e-20 LV', linestyle='--') + +print '1e20 LV min_llh', golem.MinLLH().likelihood +print '1e20 LV expectation', exp + +npp = gf.NewPhysicsParams() +npp.type = gf.NewPhysicsType.LorentzViolation +npp.n_lv = 1 +npp.lambda_1 = 1.e10 +npp.lambda_2 = 1.e10 + +golem.SetNewPhysicsParams(npp) + +exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, + drawstyle='steps-pre', label='1e-10 LV', linestyle='--') + +print '1e10 LV min_llh', golem.MinLLH().likelihood +print '1e10 LV expectation', exp + +npp = gf.NewPhysicsParams() +npp.type = gf.NewPhysicsType.LorentzViolation +npp.n_lv = 1 +npp.lambda_1 = 1.e-20 +npp.lambda_2 = 1.e-20 + +golem.SetNewPhysicsParams(npp) + +ax.tick_params(axis='x', labelsize=12) +ax.tick_params(axis='y', labelsize=12) +ax.set_xlabel(r'Deposited energy / GeV') +ax.set_ylabel(r'Events') +for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1) +for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1) + +legend = ax.legend(prop=dict(size=12)) +fig.savefig('test.png', bbox_inches='tight', dpi=250) + diff --git a/test/test_NSI.py b/test/test_NSI.py new file mode 100644 index 0000000..d144420 --- /dev/null +++ b/test/test_NSI.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python + +from __future__ import absolute_import, division + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rc + +import GolemFitPy as gf + +rc('text', usetex=True) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +dp = gf.DataPaths() +steer = gf.SteeringParams() +npp = gf.NewPhysicsParams() + +steer.quiet = False +steer.fastmode = True + +golem = gf.GolemFit(dp, steer, npp) + +fit_params = gf.FitParameters(gf.sampleTag.HESE) +golem.SetupAsimov(fit_params) + +fig = plt.figure(figsize=[6, 5]) +ax = fig.add_subplot(111) +ax.set_xscale('log') +ax.set_yscale('log') + +binning = golem.GetEnergyBinsMC() +ax.set_xlim(binning[0], binning[-1]) +# ax.set_ylim(binning[0], binning[-1]) + +print 'NULL min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='NULL', linestyle='--') + +# print 'NULL expectation', exp +print + +npp.type = gf.NewPhysicsType.NonStandardInteraction +npp.epsilon_mutau = 0.1 +golem.SetNewPhysicsParams(npp) + +print '0.1 mutau min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='0.1 mutau', linestyle='--') + +# print '0.1 mutau expectation', exp +print + +np.epsilon_mutau = 0.2 +golem.SetNewPhysicsParams(npp) + +print '0.2 mutau min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='0.2 mutau', linestyle='--') + +# print '0.2 mutau expectation', exp +print + +np.epsilon_mutau = 0.3 +golem.SetNewPhysicsParams(npp) + +print '0.3 mutau min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='0.3 mutau', linestyle='--') + +# print '0.3 mutau expectation', exp +print + +np.epsilon_mutau = 0.4 +golem.SetNewPhysicsParams(npp) + +print '0.4 mutau min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='0.4 mutau', linestyle='--') + +# print '0.4 mutau expectation', exp +print + +ax.tick_params(axis='x', labelsize=12) +ax.tick_params(axis='y', labelsize=12) +ax.set_xlabel(r'Deposited energy / GeV') +ax.set_ylabel(r'Events') +for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1) +for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1) + +legend = ax.legend(prop=dict(size=12)) +fig.savefig('test_NSI.png', bbox_inches='tight', dpi=250) + diff --git a/test/test_all_gf.py b/test/test_all_gf.py new file mode 100644 index 0000000..04b0980 --- /dev/null +++ b/test/test_all_gf.py @@ -0,0 +1,79 @@ +import numpy as np +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt + +import GolemFitPy as gf + +FASTMODE = True +PARAMETERS = [ + # 'astroFlavorAngle1', 'astroFlavorAngle2', + 'convNorm', + # 'promptNorm', 'muonNorm', 'astroNorm' +] +DEFAULTS = [ + # 4/9., 0., 1., 0., 1., 6.9 + 1. +] +RANGES = [ + # (0, 1), (-1, 1), (0.01, 10), (0., 30), (0.01, 10), (0.01, 30) + (0.01, 10) +] +BINS = 50 + +def steering_params(): + steering_categ = 'p2_0' + params = gf.SteeringParams(gf.sampleTag.HESE) + if FASTMODE: + params.fastmode = True + params.quiet = False + params.simToLoad= steering_categ.lower() + return params + +def set_up_as(fitter, params): + print 'Injecting the model', params + asimov_params = gf.FitParameters(gf.sampleTag.HESE) + for x in params.iterkeys(): + asimov_params.__setattr__(x, float(params[x])) + fitter.SetupAsimov(asimov_params) + priors = gf.Priors() + priors.convNormWidth = 9e9 + fitter.SetFitPriors(priors) + +def setup_fitter(asimov_paramset): + datapaths = gf.DataPaths() + sparams = steering_params() + npp = gf.NewPhysicsParams() + fitter = gf.GolemFit(datapaths, sparams, npp) + set_up_as(fitter, asimov_paramset) + return fitter + +def get_llh(fitter, params): + fitparams = gf.FitParameters(gf.sampleTag.HESE) + for x in params.iterkeys(): + fitparams.__setattr__(x, float(params[x])) + llh = -fitter.EvalLLH(fitparams) + return llh + +for ip, param in enumerate(PARAMETERS): + asimov_paramset = {param: DEFAULTS[ip]} + print 'injecting', asimov_paramset + fitter = setup_fitter(asimov_paramset) + binning = np.linspace(RANGES[ip][0], RANGES[ip][1], BINS) + llhs = [] + for b in binning: + test_paramset = {param: b} + print 'testing', test_paramset + llh = get_llh(fitter, test_paramset) + print 'llh', llh + llhs.append(llh) + plt.plot(binning, llhs) + plt.axvline(x=DEFAULTS[ip]) + plt.xlabel(param) + plt.ylabel('LLH') + outfile = 'llh_profile_noprior_' + if FASTMODE: + plt.savefig(outfile + 'fastmode_{0}.png'.format(param)) + else: + plt.savefig(outfile + '{0}.png'.format(param)) + plt.clf() diff --git a/test/test_gf.py b/test/test_gf.py index 6f8f463..4b77924 100644 --- a/test/test_gf.py +++ b/test/test_gf.py @@ -1,9 +1,14 @@ +import numpy.ma as ma + import GolemFitPy as gf +FASTMODE = True + def steering_params(): steering_categ = 'p2_0' params = gf.SteeringParams(gf.sampleTag.HESE) - params.fastmode = True + if FASTMODE: + params.fastmode = True params.quiet = False params.simToLoad= steering_categ.lower() return params @@ -50,7 +55,7 @@ test_paramset = {'astroFlavorAngle1': 4/9., 'astroFlavorAngle2': 0} print 'testing', test_paramset print 'llh', get_llh(fitter, test_paramset) -test_paramset = {'astroFlavorAngle1': 4/9, 'astroFlavorAngle2': 0} +test_paramset = {'astroFlavorAngle1': 4/9., 'astroFlavorAngle2': 0} print 'testing', test_paramset print 'llh', get_llh(fitter, test_paramset) @@ -68,9 +73,17 @@ for an1 in angle1_binning: llhs = [] for an2 in angle2_binning: test_paramset = {'astroFlavorAngle1': an1, 'astroFlavorAngle2': an2} - llhs.append(get_llh(fitter, test_paramset)) + llh = get_llh(fitter, test_paramset) + if abs(llh) > 9e9: + llhs.append(np.nan) + else: + llhs.append(llh) + llhs = ma.masked_invalid(llhs) plt.plot(angle2_binning, llhs, label='astroFlavorAngle1 = {0}'.format(an1)) plt.xlabel('astroFlavorAngle2') plt.ylabel('LLH') plt.legend() -plt.savefig('llh_profile_fastmode.png'.format(an1)) +if FASTMODE: + plt.savefig('llh_profile_fastmode.png'.format(an1)) +else: + plt.savefig('llh_profile.png'.format(an1)) diff --git a/test/test_gf_simple.py b/test/test_gf_simple.py index 51460c1..3a6ebb5 100644 --- a/test/test_gf_simple.py +++ b/test/test_gf_simple.py @@ -1,11 +1,20 @@ +import numpy as np +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt + import GolemFitPy as gf +FASTMODE = True + dp = gf.DataPaths() npp = gf.NewPhysicsParams() sp = gf.SteeringParams(gf.sampleTag.HESE) sp.quiet = False -# sp.fastmode = True +if FASTMODE: + sp.fastmode = True +sp.frequentist = True golem = gf.GolemFit(dp, sp, npp) @@ -20,4 +29,33 @@ fp_sh.astroFlavorAngle1 = 0.36 fp_sh.astroFlavorAngle2 = -0.57 print 'Eval fp = {0}'.format(golem.EvalLLH(fp)) + +# energy_centers = golem.GetEnergyBinsMC()[:-1]+ np.diff(golem.GetEnergyBinsMC())/2. + +# plt.hist(energy_centers,bins=golem.GetEnergyBinsMC(), +# weights=np.sum(golem.GetExpectation(fp),axis=(0,1,2,3)), +# histtype="step", lw = 2, label='injected') + +# data_energy_dist = np.sum(golem.GetDataDistribution(),axis=(0,1,2,3)) +# energy_centers=golem.GetEnergyBinsData()[:-1]+ np.diff(golem.GetEnergyBinsData())/2. +# plt.errorbar(energy_centers,data_energy_dist,yerr = np.sqrt(data_energy_dist),fmt='o') + print 'Eval fp_sh = {0}'.format(golem.EvalLLH(fp_sh)) + +# plt.hist(energy_centers,bins=golem.GetEnergyBinsMC(), +# weights=np.sum(golem.GetExpectation(fp_sh),axis=(0,1,2,3)), +# histtype="step", lw = 2, label='test') + +# data_energy_dist = np.sum(golem.GetDataDistribution(),axis=(0,1,2,3)) +# energy_centers=golem.GetEnergyBinsData()[:-1]+ np.diff(golem.GetEnergyBinsData())/2. +# plt.errorbar(energy_centers,data_energy_dist,yerr = np.sqrt(data_energy_dist),fmt='o') + +# plt.loglog(nonposy="clip") +# plt.xlabel(r"Deposited energy/GeV") +# plt.ylabel(r"Events") + +# outname = 'Expectation' +# if FASTMODE: +# plt.savefig(outname + 'fastmode.png') +# else: +# plt.savefig(outname + '.png') diff --git a/utils/fr.py b/utils/fr.py index 61b5a61..5acd6f8 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -212,7 +212,7 @@ def estimate_scale(args): upper_s = (m_eign/args.binning[0]) / (args.binning[0]**(args.dimension-3)) scale = np.power(10, np.average(np.log10([lower_s, upper_s]))) diff = upper_s / lower_s - scale_region = (lower_s/diff, upper_s*diff) + scale_region = (lower_s/np.power(10, args.dimension), upper_s*diff*np.power(10, args.dimension)) scale_region = [np.power(10, np.round(np.log10(x))) for x in scale_region] return scale, scale_region diff --git a/utils/gf.py b/utils/gf.py index ea4f61f..098faa9 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -57,9 +57,14 @@ def steering_params(args): steering_categ = args.ast params = gf.SteeringParams() params.quiet = False + params.fastmode = True params.simToLoad= steering_categ.name.lower() params.evalThreads = args.threads # params.evalThreads = thread_factors(args.threads)[1] + + # For Tianlu + # params.years = [999] + return params diff --git a/utils/likelihood.py b/utils/likelihood.py index 2e4a22d..c2822a1 100644 --- a/utils/likelihood.py +++ b/utils/likelihood.py @@ -167,9 +167,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): param.value = flavour_angles[idx] - print 'llh_paramset', llh_paramset - print 'hypo_paramset', hypo_paramset - print 'fr', fr if args.likelihood is Likelihood.FLAT: llh = 1. elif args.likelihood is Likelihood.GAUSSIAN: @@ -179,7 +176,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): llh = gf_utils.get_llh(fitter, hypo_paramset) elif args.likelihood is Likelihood.GF_FREQ: lhh = gf_utils.get_llh_freq(fitter, hypo_paramset) - print 'llh', llh return llh diff --git a/utils/misc.py b/utils/misc.py index cad03bc..5483675 100644 --- a/utils/misc.py +++ b/utils/misc.py @@ -12,6 +12,7 @@ from __future__ import absolute_import, division import os import errno import multiprocessing +from fractions import gcd import argparse from operator import attrgetter @@ -33,14 +34,18 @@ class SortingHelpFormatter(argparse.HelpFormatter): super(SortingHelpFormatter, self).add_arguments(actions) +def solve_ratio(fr): + denominator = reduce(gcd, fr) + return [int(x/denominator) for x in fr] + + def gen_identifier(args): f = '_DIM{0}'.format(args.dimension) - mr1, mr2, mr3 = args.measured_ratio + mr1, mr2, mr3 = solve_ratio(args.measured_ratio) if args.fix_source_ratio: - sr1, sr2, sr3 = args.source_ratio - f += '_sfr_{0:03d}_{1:03d}_{2:03d}_mfr_{3:03d}_{4:03d}_{5:03d}'.format( - int(sr1*100), int(sr2*100), int(sr3*100), - int(mr1*100), int(mr2*100), int(mr3*100) + sr1, sr2, sr3 = solve_ratio(args.source_ratio) + f += '_sfr_{0:G}_{1:G}_{2:G}_mfr_{3:G}_{4:G}_{5:G}'.format( + sr1, sr2, sr3, mr1, mr2, mr3 ) if args.fix_mixing: f += '_fix_mixing' @@ -128,6 +133,21 @@ def make_dir(outfile): else: raise +def remove_option(parser, arg): + for action in parser._actions: + if (vars(action)['option_strings'] + and vars(action)['option_strings'][0] == arg) \ + or vars(action)['dest'] == arg: + parser._remove_action(action) + + for action in parser._action_groups: + vars_action = vars(action) + var_group_actions = vars_action['_group_actions'] + for x in var_group_actions: + if x.dest == arg: + var_group_actions.remove(x) + return + def seed_parse(s): if s.lower() == 'none': diff --git a/utils/multinest.py b/utils/multinest.py index 9dd0742..f3595f9 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 +from utils.misc import gen_identifier, make_dir def CubePrior(cube, ndim, n_params): @@ -37,13 +37,15 @@ def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset, llh_paramset[pm].value = mn_paramset[pm].value theta = llh_paramset.values # print 'llh_paramset', llh_paramset - return likelihood.ln_prob( + llh = likelihood.ln_prob( theta=theta, args=args, asimov_paramset=asimov_paramset, llh_paramset=llh_paramset, fitter=fitter ) + # print 'llh', llh + return llh def mn_argparse(parser): @@ -77,8 +79,8 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter): fitter = fitter ) - prefix = './mnrun/DIM{0}/{1}_{2:>019}_'.format( - args.dimension, gen_outfile_name(args), np.random.randint(0, 2**63) + prefix = './mnrun/DIM{0}/{1}_{2:>010}_'.format( + args.dimension, gen_identifier(args), np.random.randint(0, 2**30) ) make_dir(prefix) print 'Running evidence calculation for {0}'.format(prefix) diff --git a/utils/param.py b/utils/param.py index f861264..25558c0 100644 --- a/utils/param.py +++ b/utils/param.py @@ -252,10 +252,16 @@ def get_paramsets(args, nuisance_paramset): ]) if not args.fix_scale: tag = ParamTag.SCALE - llh_paramset.append( - Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag) - ) + if hasattr(args, 'dimension'): + llh_paramset.append( + Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, + tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag) + ) + elif hasattr(args, 'dimensions'): + llh_paramset.append( + Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, + tex=r'{\rm log}_{10}\Lambda^{-1} / GeV^{-d+4}', tag=tag) + ) if not args.fix_source_ratio: tag = ParamTag.SRCANGLES llh_paramset.extend([ diff --git a/utils/plot.py b/utils/plot.py index 3d94cc1..8f4eba8 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -10,6 +10,7 @@ Plotting functions for the BSM flavour ratio analysis from __future__ import absolute_import, division import os +from copy import deepcopy import numpy as np import matplotlib as mpl @@ -112,18 +113,18 @@ def plot_Tchain(Tchain, axes_labels, ranges): def gen_figtext(args): """Generate the figure text.""" t = '' - mr1, mr2, mr3 = args.measured_ratio + mr1, mr2, mr3 = misc_utils.solve_ratio(args.measured_ratio) if args.fix_source_ratio: - sr1, sr2, sr3 = args.source_ratio - t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \ - 'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \ - '{5:.2f}]\nDimension = {6}'.format( + 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) ) else: - t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ - '{2:.2f}]\nDimension = {3}'.format( + t += 'IC observed flavour ratio = [{0}, {1}, ' \ + '{2}]\nDimension = {3}'.format( mr1, mr2, mr3, args.dimension, int(args.energy) ) if args.fix_scale: @@ -229,28 +230,33 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset): g.export(outfile+'_elements.'+of) +def myround(x, base=5, 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)) + else: int(base * np.round(float(x)/base)) + + def plot_statistic(data, outfile, outformat, args, scale_param, label=None): """Make MultiNest factor or LLH value plot.""" - print "Making Statistic plot" + print 'Making Statistic plot' fig_text = gen_figtext(args) if label is not None: fig_text += '\n' + label print 'data', data print 'data.shape', data.shape scales, statistic = data.T + min_idx = np.argmin(scales) + null = statistic[min_idx] if args.stat_method is StatCateg.BAYESIAN: - min_idx = np.argmin(scales) - null = statistic[min_idx] reduced_ev = -(statistic - null) elif args.stat_method is StatCateg.FREQUENTIST: - min_idx = np.argmin(scales) - null = statistic[min_idx] reduced_ev = -2*(statistic - null) fig = plt.figure(figsize=(7, 5)) ax = fig.add_subplot(111) - ax.set_xlim(scale_param.ranges) + ax.set_xlim(np.log10(args.scale_region)) ax.set_xlabel('$'+scale_param.tex+'$') if args.stat_method is StatCateg.BAYESIAN: ax.set_ylabel(r'Bayes Factor') @@ -271,85 +277,302 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None): ax.add_artist(at) misc_utils.make_dir(outfile) for of in outformat: + print 'Saving as {0}'.format(outfile+'.'+of) fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) -def myround(x, base=5, 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)) - else: int(base * np.round(float(x)/base)) +def plot_sens_full(data, outfile, outformat, args): + print 'Making FULL sensitivity plot' - -def plot_BSM_angles_limit(dirname, outfile, outformat, args, bayesian): - """Make BSM angles vs scale limit plot.""" - if not args.plot_angles_limit: return - print "Making BSM angles limit plot.""" - fig_text = gen_figtext(args) - xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$'] - - raw = [] - for root, dirs, filenames in os.walk(dirname): - for fn in filenames: - if fn[-4:] == '.npy': - raw.append(np.load(os.path.join(root, fn))) - raw = np.vstack(raw) - raw = raw[np.argsort(raw[:,0])] - print 'raw', raw - print 'raw.shape', raw.shape - sc_ranges = ( - myround(np.min(raw[0][:,0]), up=True), - myround(np.max(raw[0][:,0]), down=True) - ) - - proc = [] - if bayesian: - for idx, theta in enumerate(raw): - scale, evidences = theta.T - null = evidences[0] - reduced_ev = -(evidences - null) - al = scale[reduced_ev > np.log(10**(1/2.))] - proc.append((idx+1, al[0])) - else: - for idx, theta in enumerate(raw): - scale, llh = theta.T - delta_llh = -2 * (llh - np.max(llh)) - # 90% CL for 1 dof - al = scale[delta_llh > 2.71] - proc.append((idx+1, al[0])) - - limits = np.array(proc) - print 'limits', limits + colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + xticks = ['{0}'.format(x) for x in range(np.min(args.dimensions), + np.max(args.dimensions)+1)] + yranges = [np.inf, -np.inf] + legend_handles = [] fig = plt.figure(figsize=(7, 5)) ax = fig.add_subplot(111) - ax.set_xticklabels(['']+xticks+['']) - ax.set_xlim(0, len(xticks)+1) - ax.set_ylim(sc_ranges[0], sc_ranges[-1]) - ax.set_xlabel(r'BSM angle') - ylabel = r'${\rm log}_{10} \Lambda' + get_units(args.dimension) + r'$' - ax.set_ylabel(ylabel) + ax.set_xlim(args.dimensions[0]-1, args.dimensions[-1]+1) + ax.set_xticklabels([''] + xticks + ['']) + ax.set_xlabel(r'BSM operator dimension ' + r'$d$') + ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$') + + argsc = deepcopy(args) + for idim in xrange(len(data)): + dim = args.dimensions[idim] + print 'dim', dim + argsc.dimension = dim + for isrc in xrange(len(data[idim])): + src = args.source_ratios[isrc] + argsc.source_ratio = src + + # fig_text = gen_figtext(argsc) + scales, statistic = data[idim][isrc].T + min_idx = np.argmin(scales) + null = statistic[min_idx] + if args.stat_method is StatCateg.BAYESIAN: + reduced_ev = -(statistic - null) + # al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief + al = scales[reduced_ev > 0.4] # Testing + elif args.stat_method is StatCateg.FREQUENTIST: + reduced_ev = -2*(statistic - null) + al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks + if len(al) == 0: + print 'No points for DIM {0} FRS {1} NULL {2}!'.format( + dim, misc_utils.solve_ratio(src), null + ) + print 'Reduced EV {0}'.format(reduced_ev) + continue + lim = al[0] + label = '[{0}, {1}, {2}]'.format(*misc_utils.solve_ratio(src)) + if lim < yranges[0]: yranges[0] = lim + if lim > yranges[1]: yranges[1] = lim+4 + line = plt.Line2D( + (dim-0.1, dim+0.1), (lim, lim), lw=3, color=colour[isrc], label=label + ) + ax.add_line(line) + if idim == 0: legend_handles.append(line) + x_offset = isrc*0.05 - 0.05 + ax.annotate( + s='', xy=(dim+x_offset, lim), xytext=(dim+x_offset, lim+3), + arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[isrc]} + ) - for l in limits: - line = plt.Line2D( - (l[0]-0.1, l[0]+0.1), (l[1], l[1]), lw=3, color='r' - ) - ax.add_line(line) - # ax.arrow( - # l[0], l[1], 0, -1.5, head_width=0.05, head_length=0.2, fc='r', - # ec='r', lw=2 - # ) - ax.annotate( - s='', xy=l, xytext=(l[0], l[1]+1.5), - arrowprops={'arrowstyle': '<-', 'lw': 1.5, 'color':'r'} - ) + try: + yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) + ax.set_ylim(yranges) + except: pass + ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right') for ymaj in ax.yaxis.get_majorticklocs(): - ax.axhline(y=ymaj, ls='-', color='gray', alpha=0.4, linewidth=1) + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) for xmaj in ax.xaxis.get_majorticklocs(): - ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1) + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + misc_utils.make_dir(outfile) for of in outformat: + print 'Saving plot as {0}'.format(outfile+'.'+of) fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) + +def plot_sens_fixed_angle(data, outfile, outformat, args): + print 'Making FIXED_ANGLE sensitivity plot' + + colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$'] + argsc = deepcopy(args) + for idim in xrange(len(data)): + dim = args.dimensions[idim] + print '= dim', dim + argsc.dimension = dim + + yranges = [np.inf, -np.inf] + legend_handles = [] + + 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} \Lambda^{-1} / GeV^{-d+4}$') + + for isrc in xrange(len(data[idim])): + src = args.source_ratios[isrc] + argsc.source_ratio = src + print '== src', src + + for ian in xrange(len(data[idim][isrc])): + print '=== an', ian + scales, statistic = data[idim][isrc][ian].T + min_idx = np.argmin(scales) + null = statistic[min_idx] + if args.stat_method is StatCateg.BAYESIAN: + reduced_ev = -(statistic - null) + al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief + elif args.stat_method is StatCateg.FREQUENTIST: + reduced_ev = -2*(statistic - null) + al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks + print 'reduced_ev', reduced_ev + if len(al) == 0: + print 'No points for DIM {0} FRS {1}\nreduced_ev {2}!'.format( + dim, src, reduced_ev + ) + continue + 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[1]: yranges[1] = lim+4 + line = plt.Line2D( + (ian+1-0.1, ian+1+0.1), (lim, lim), lw=3, color=colour[isrc], label=label + ) + ax.add_line(line) + if len(legend_handles) < isrc+1: + legend_handles.append(line) + x_offset = isrc*0.05 - 0.05 + ax.annotate( + s='', xy=(ian+1+x_offset, lim), xytext=(ian+1+x_offset, lim+3), + arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[isrc]} + ) + + try: + yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) + ax.set_ylim(yranges) + except: pass + + ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right') + for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) + for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + + out = outfile + '_DIM{0}'.format(dim) + misc_utils.make_dir(out) + for of in outformat: + print 'Saving plot as {0}'.format(out+'.'+of) + fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150) + + + for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) + for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + + +def plot_sens_corr_angle(data, outfile, outformat, args): + print 'Making CORR_ANGLE sensitivity plot' + + labels = [r'$sin^2(2\mathcal{O}_{12})$', + r'$sin^2(2\mathcal{O}_{13})$', + r'$sin^2(2\mathcal{O}_{23})$'] + + argsc = deepcopy(args) + for idim in xrange(len(data)): + dim = args.dimensions[idim] + print '= dim', dim + argsc.dimension = dim + for isrc in xrange(len(data[idim])): + src = args.source_ratios[isrc] + argsc.source_ratio = src + print '== src', src + for ian in xrange(len(data[idim][isrc])): + print '=== an', ian + + d = data[idim][isrc][ian].reshape( + len(data[idim][isrc][ian])**2, 3 + ) + + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + ax.set_ylim(0, 1) + ax.set_ylabel(labels[ian]) + ax.set_xlabel(r'${\rm log}_{10} \Lambda^{-1}'+get_units(dim)+r'$') + + xranges = [np.inf, -np.inf] + legend_handles = [] + + y = d[:,0] + x = d[:,1] + z = d[:,2] + + print 'x', x + print 'y', y + print 'z', z + null_idx = np.argmin(x) + null = z[null_idx] + print 'null = {0}, for scale = {1}'.format(null, x[null_idx]) + + if args.stat_method is StatCateg.BAYESIAN: + z_r = -(z - null) + elif args.stat_method is StatCateg.FREQUENTIST: + z_r = -2*(z - null) + print 'scale', x + print 'reduced ev', z_r + + pdat = np.array([x, y, z_r, np.ones(x.shape)]).T + print 'pdat', pdat + pdat_clean = [] + for d in pdat: + if not np.any(np.isnan(d)): pdat_clean.append(d) + pdat = np.vstack(pdat_clean) + sort_column = 3 + pdat_sorted = pdat[pdat[:,sort_column].argsort()] + uni, c = np.unique(pdat[:,sort_column], return_counts=True) + print uni, c + print len(uni) + print np.unique(c) + + n = len(uni) + assert len(np.unique(c)) == 1 + c = c[0] + col_array = [] + for col in pdat_sorted.T: + col_array.append(col.reshape(n, c)) + col_array = np.stack(col_array) + sep_arrays = [] + for x_i in xrange(n): + sep_arrays.append(col_array[:,x_i]) + + print len(sep_arrays) + sep_arrays = sep_arrays[0][:3] + print sep_arrays + + if args.stat_method is StatCateg.BAYESIAN: + reduced_pdat_mask = (sep_arrays[2] > log(10**(3/2.))) # Strong degree of belief + elif args.stat_method is StatCateg.FREQUENTIST: + reduced_pdat_mask = (sep_arrays[2] > 4.61) # 90% CL for 2 DOFS via Wilks + reduced_pdat = sep_arrays.T[reduced_pdat_mask].T + print 'reduced_pdat', reduced_pdat + + ax.tick_params(axis='x', labelsize=11) + ax.tick_params(axis='y', labelsize=11) + + mini, maxi = np.min(x), np.max(x) + ax.set_xlim((mini, maxi)) + ax.set_ylim(0, 1) + ax.grid(b=False) + + x_v = reduced_pdat[0].round(decimals=4) + y_v = reduced_pdat[1].round(decimals=4) + uniques = np.unique(x_v) + print 'uniques', uniques + if len(uniques) == 1: continue + bw = np.min(np.diff(uniques)) + print 'bw', bw + print np.diff(uniques) + uni_x_split = np.split(uniques, np.where(np.diff(uniques) > bw*(1e20))[0] + 1) + print 'len(uni_x_split)', len(uni_x_split) + for uni_x in uni_x_split: + p_x_l, p_y_l = [], [] + p_x_u, p_y_u = [], [] + for uni in uni_x: + idxes = np.where(x_v == uni)[0] + ymin, ymax = 1, 0 + for idx in idxes: + if y_v[idx] < ymin: ymin = y_v[idx] + if y_v[idx] > ymax: ymax = y_v[idx] + p_x_l.append(uni) + p_y_l.append(ymin) + p_x_u.append(uni) + p_y_u.append(ymax) + p_x_l, p_y_l = np.array(p_x_l, dtype=np.float64), np.array(p_y_l, dtype=np.float64) + p_x_u, p_y_u = np.array(list(reversed(p_x_u)), dtype=np.float64), np.array(list(reversed(p_y_u)), dtype=np.float64) + p_x = np.hstack([p_x_l, p_x_u]) + p_y = np.hstack([p_y_l, p_y_u]) + p_x = np.r_[p_x, p_x[0]] + p_y = np.r_[p_y, p_y[0]] + print 'p_x', p_x + print 'p_y', p_y + try: + tck, u = interpolate.splprep([p_x, p_y], s=1e-5, per=True) + xi, yi = interpolate.splev(np.linspace(0, 1, 1000), tck) + except: + xi, yi = p_x, p_y + ax.fill(xi, yi, 'r', edgecolor='r', linewidth=1) + + out = outfile + '_DIM{0}_SRC{1}_AN{2}'.format(dim, isrc, ian) + misc_utils.make_dir(out) + for of in outformat: + print 'Saving plot as {0}'.format(out+'.'+of) + fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150) |