diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rwxr-xr-x | fr.py | 136 | ||||
| -rw-r--r-- | mnrun/plot.py | 23 | ||||
| -rw-r--r-- | plot_bf.py | 252 | ||||
| -rwxr-xr-x | submitter/clean.sh | 2 | ||||
| -rw-r--r-- | submitter/make_dag.py | 53 | ||||
| -rw-r--r-- | submitter/submit.sub | 8 | ||||
| -rw-r--r-- | utils/plot.py | 4 |
8 files changed, 417 insertions, 62 deletions
@@ -8,3 +8,4 @@ dagman_FR* submitter/logs/ mnrun_* *.png +mnrun/ @@ -87,11 +87,9 @@ def get_paramsets(args): 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: - logLam, scale_region = np.log10(args.scale), np.log10(args.scale_region) - lL_range = (logLam-scale_region, logLam+scale_region) tag = ParamTag.SCALE mcmc_paramset.append( - Param(name='logLam', value=logLam, ranges=lL_range, std=3, + 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: @@ -140,6 +138,7 @@ def process_args(args): ) + 6 ) """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 @@ -147,7 +146,7 @@ def process_args(args): args.bayes_eval_bin = int(args.bayes_eval_bin) -def parse_args(): +def parse_args(args=None): """Parse command line arguments""" parser = argparse.ArgumentParser( description="BSM flavour ratio analysis", @@ -187,6 +186,10 @@ def parse_args(): 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' ) @@ -211,6 +214,10 @@ def parse_args(): 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' ) @@ -263,7 +270,8 @@ def parse_args(): gf_utils.gf_argparse(parser) plot_utils.plot_argparse(parser) nuisance_argparse(parser) - return parser.parse_args() + if args is None: return parser.parse_args() + else: return parser.parse_args(args.split()) def main(): @@ -320,9 +328,8 @@ def main(): ) out = args.bayes_output+'/fr_evidence' + misc_utils.gen_identifier(args) - sc_range = mcmc_paramset.from_tag(ParamTag.SCALE)[0].ranges scan_scales = np.linspace( - sc_range[0], sc_range[1], args.bayes_bins + np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.bayes_bins ) if args.run_bayes_factor: import pymultinest @@ -362,8 +369,9 @@ def main(): fitter=fitter ) - prefix = 'mnrun_{0:.0E}'.format(np.power(10, sc)) + \ - '_' + os.path.basename(out) + '_' + prefix = 'mnrun/' + os.path.basename(out) + '_' + if args.bayes_eval_bin is not None: + prefix += '{0}_'.format(s_idx) print 'begin running evidence calculation for {0}'.format(prefix) result = pymultinest.run( LogLikelihood=lnProb, @@ -389,8 +397,7 @@ def main(): dirname = os.path.dirname(out) plot_utils.bayes_factor_plot( - dirname=dirname, outfile=out, outformat=['png'], args=args, - xlim=sc_range + dirname=dirname, outfile=out, outformat=['png'], args=args ) out = args.angles_lim_output+'/fr_an_evidence' + misc_utils.gen_identifier(args) @@ -414,7 +421,9 @@ def main(): # default are uniform priors return ; - data = np.zeros((len(scenarios), args.bayes_bins, 2)) + if args.bayes_eval_bin is not None: + data = np.zeros((len(scenarios), 1, 2)) + else: data = np.zeros((len(scenarios), args.bayes_bins, 2)) mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES) sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0] for idx, scen in enumerate(scenarios): @@ -446,8 +455,9 @@ def main(): mcmc_paramset=mcmc_paramset, fitter=fitter ) - prefix = 'mnrun_{0:.0E}'.format(np.power(10, sc)) + \ - '_' + os.path.basename(out) + '_' + prefix = 'mnrun/' + os.path.basename(out) + '_' + if args.bayes_eval_bin is not None: + prefix += '{0}_{1}_'.format(idx, s_idx) print 'begin running evidence calculation for {0}'.format(prefix) result = pymultinest.run( LogLikelihood=lnProb, @@ -481,6 +491,104 @@ def main(): args=args, bayesian=True ) + out = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args) + if args.run_angles_correlation: + if args.bayes_eval_bin is None: assert 0 + import pymultinest + + scenarios = [ + [1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 1, 0], + ] + nuisance = mcmc_paramset.from_tag(ParamTag.NUISANCE) + mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES) + sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0] + + if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT: + fitter = gf_utils.setup_fitter(args, asimov_paramset) + else: fitter = None + + def CubePrior(cube, ndim, nparams): + # default are uniform priors + return ; + + scan_angles = np.linspace(0, 1, args.bayes_bins) + + if args.bayes_eval_bin is not None: + data = np.zeros((len(scenarios), 1, 1, 3)) + else: data = np.zeros((len(scenarios), scale_bins, angle_bins, 3)) + for idx, scen in enumerate(scenarios): + for an in mm_angles: + an.value = 0 + keep = mcmc_paramset.from_tag(ParamTag.MMANGLES)[idx] + q = ParamSet(nuisance.params + [x for x in mm_angles if x.name != keep.name]) + n_params = len(q) + prior_ranges = q.seeds + + scales, angles, evidences = [], [], [] + for s_idx, sc in enumerate(scan_scales): + for a_idx, an in enumerate(scan_angles): + index = s_idx*args.bayes_bins + a_idx + if args.bayes_eval_bin is not None: + if args.bayes_eval_bin == index: + if idx == 0: + out += '_scale_{0:.0E}'.format(np.power(10, sc)) + else: continue + + print '== SCALE = {0:.0E}'.format(np.power(10, sc)) + print '== ANGLE = {0:.2f}'.format(an) + sc_angles.value = sc + keep.value = an + def lnProb(cube, ndim, nparams): + for i in range(ndim-1): + prange = prior_ranges[i][1] - prior_ranges[i][0] + q[i].value = prange*cube[i] + prior_ranges[i][0] + for name in q.names: + mcmc_paramset[name].value = q[name].value + theta = mcmc_paramset.values + # print 'theta', theta + # print 'mcmc_paramset', mcmc_paramset + return llh_utils.triangle_llh( + theta=theta, + args=args, + asimov_paramset=asimov_paramset, + mcmc_paramset=mcmc_paramset, + fitter=fitter + ) + prefix = 'mnrun/' + os.path.basename(out) + '_' + if args.bayes_eval_bin is not None: + prefix += '{0}_{1}_{2}'.format(idx, s_idx, a_idx) + + print 'begin running evidence calculation for {0}'.format(prefix) + result = pymultinest.run( + LogLikelihood=lnProb, + Prior=CubePrior, + n_dims=n_params, + importance_nested_sampling=True, + n_live_points=args.bayes_live_points, + evidence_tolerance=args.bayes_tolerance, + outputfiles_basename=prefix, + resume=False, + verbose=True + ) + + analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params) + a_lnZ = analyzer.get_stats()['global evidence'] + print 'Evidence = {0}'.format(a_lnZ) + scales.append(sc) + angles.append(an) + evidences.append(a_lnZ) + + for i, d in enumerate(evidences): + data[idx][i][i][0] = scales[i] + data[idx][i][i][1] = angles[i] + data[idx][i][i][2] = d + + misc_utils.make_dir(out) + print 'saving to {0}.npy'.format(out) + np.save(out+'.npy', np.array(data)) + print "DONE!" diff --git a/mnrun/plot.py b/mnrun/plot.py deleted file mode 100644 index f7e750b..0000000 --- a/mnrun/plot.py +++ /dev/null @@ -1,23 +0,0 @@ -import numpy as np -import matplotlib as mpl -mpl.use('Agg') -from matplotlib import pyplot as plt -from matplotlib import rc - -rc('text', usetex=False) -rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) - -arr = np.load('fr_evidence_test_033_033_033_0010_sfr_033_066_000_DIM6_single_scale.npy') - -print arr -ev = arr.T[1] -null = ev[0] -print null -re_z = -(ev - null) -print re_z - -plt.plot(arr.T[0], re_z) -plt.xlabel(r'${\rm log}_{10} \Lambda$') -plt.ylabel(r'Bayes Factor') -# plt.axhline(np.power(10, 1.5), color='grey', linestyle='-', alpha=0.4) -plt.savefig('./test.png', bbox_inches='tight', dpi=150) diff --git a/plot_bf.py b/plot_bf.py new file mode 100644 index 0000000..faefcde --- /dev/null +++ b/plot_bf.py @@ -0,0 +1,252 @@ +#! /usr/bin/env python +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : April 14, 2018 + +""" +HESE BSM flavour ratio sensivity plotting script +""" + +from __future__ import absolute_import, division + +import os + +import numpy as np +import numpy.ma as ma + +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import rc +from matplotlib import pyplot as plt +from matplotlib.offsetbox import AnchoredText + +import fr +from utils import misc as misc_utils +from utils.fr import normalise_fr +from utils.plot import myround, get_units + + +rc('text', usetex=False) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +fix_sfr_mfr = [ + (1, 1, 1, 1, 2, 0), + # (1, 1, 1, 1, 0, 0), + # (1, 1, 1, 0, 1, 0), +] + +# FR +dimension = [3, 4, 5, 6, 7, 8] +sigma_ratio = ['0.01'] +energy_dependance = 'spectral' +spectral_index = -2 +binning = [1e4, 1e7, 5] +fix_mixing = 'False' +fix_mixing_almost = 'False' +scale_region = "1E10" + +# Likelihood +likelihood = 'golemfit' + +# Nuisance +convNorm = 1. +promptNorm = 0. +muonNorm = 1. +astroNorm = 6.9 +astroDeltaGamma = 2.5 + +# GolemFit +ast = 'p2_0' +data = 'real' + +# Bayes Factor +bayes_bins = 100 +bayes_live_points = 1000 +bayes_tolerance = 0.01 +bayes_eval_bin = True # set to 'all' to run normally + +# Plot +plot_bayes = False +plot_angles_limit = False +plot_angles_corr = True +outformat = ['png'] +significance = np.log(10**(1/2.)) +# significance = np.log(10**(3/2.)) + + +bayes_array = ma.masked_equal(np.zeros((len(dimension), len(fix_sfr_mfr), bayes_bins, 2)), 0) +angles_lim_array = np.zeros((len(dimension), len(fix_sfr_mfr), 3, bayes_bins, 2)) +for i_dim, dim in enumerate(dimension): + if energy_dependance == 'mono': + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/{2:.0E}'.format(likelihood, dim, en) + elif energy_dependance == 'spectral': + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(likelihood, dim, spectral_index) + + bayes_output = 'None' + angles_lim_output = 'None' + angles_corr_output = 'None' + for sig in sigma_ratio: + for i_frs, frs in enumerate(fix_sfr_mfr): + outchains = outchain_head + '/fix_ifr/{0}/'.format(str(sig).replace('.', '_')) + if plot_bayes: + bayes_output = outchains + '/bayes_factor/' + if plot_angles_limit: + angles_lim_output = outchains + '/angles_limit/' + if plot_angles_corr: + angles_corr_output = outchains + '/angles_corr/' + + argstring = '--measured-ratio {0} {1} {2} --fix-source-ratio True --source-ratio {3} {4} {5} --dimension {6} --seed 24 --outfile {7} --run-mcmc False --likelihood {8} --plot-angles False --bayes-output {9} --angles-lim-output {10} --bayes-bins {11} --angles-corr-output'.format(frs[0], frs[1], frs[2], frs[3], frs[4], frs[5], dim, outchains, likelihood, bayes_output, angles_lim_output, bayes_bins, angles_corr_output) + args = fr.parse_args(argstring) + fr.process_args(args) + # misc_utils.print_args(args) + + if plot_bayes: + infile = args.bayes_output+'/fr_evidence'+misc_utils.gen_identifier(args) + if plot_angles_limit: + infile = args.angles_lim_output+'/fr_an_evidence'+misc_utils.gen_identifier(args) + if plot_angles_corr: + infile = args.angles_corr_output+'/fr_co_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 + ) + raw = [] + fail = 0 + for sc in scan_scales: + try: + infile_s = infile + '_scale_{0:.0E}'.format(np.power(10, sc)) + lf = np.load(infile_s+'.npy') + if plot_angles_limit: + if len(lf.shape) == 3: lf = lf[:,0,:] + if plot_angles_corr: + # TODO(shivesh) + assert 0 + if len(lf.shape) == 3: lf = lf[:,0,:] + raw.append(lf) + except IOError: + fail += 1 + print 'failed to open {0}'.format(infile_s) + if plot_bayes: + raw.append([0, 0]) + if plot_angles_limit: + raw.append(np.zeros((3, 2))) + pass + print 'failed to open {0} files'.format(fail) + + if plot_bayes: + raw = np.vstack(raw) + if plot_angles_limit: + raw = np.vstack(raw).reshape(args.bayes_bins, 3, 2) + a = ma.masked_equal(np.zeros((3, args.bayes_bins, 2)), 0) + for i_x, x in enumerate(raw): + for i_y, y in enumerate(x): + a[i_y][i_x] = ma.masked_equal(y, 0) + + if plot_bayes: + bayes_array[i_dim][i_frs] = ma.masked_equal(raw, 0) + + if plot_angles_limit: + # TODO(shivesh) + angles_lim_array[i_dim][i_frs] = ma.masked_equal(a, 0) + +if plot_bayes: + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + + colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + yranges = [np.inf, -np.inf] + legend_handles = [] + ax.set_xlim(dimension[0]-1, dimension[-1]+1) + xticks = [''] + range(dimension[0], dimension[-1]+1) + [''] + ax.set_xticklabels(xticks) + ax.set_xlabel(r'BSM operator dimension ' + r'$d$') + ax.set_ylabel(r'${\rm log}_{10} \Lambda / GeV^{-d+4}$') + for i_dim, dim in enumerate(dimension): + for i_frs, frs in enumerate(fix_sfr_mfr): + scale, evidences = bayes_array[i_dim][i_frs].T + null = evidences[np.argmin(scale)] + # TODO(shivesh): negative or not? + reduced_ev = -(evidences - null) + al = scale[reduced_ev > significance] + if len(al) > 0: + label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5]) + lim = al[0] + print 'frs, dim, lim = ', frs, dim, lim + 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[i_frs], label=label + ) + ax.add_line(line) + if i_dim == 0: legend_handles.append(line) + x_offset = i_frs*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[i_frs]} + ) + + else: + print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, null) + yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True)) + ax.set_ylim(yranges) + + 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) + + for of in outformat: + fig.savefig('./images/bayes_factor.'+of, bbox_inches='tight', dpi=150) + +if plot_angles_limit: + colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} + for i_dim, dim in enumerate(dimension): + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + yranges = [np.inf, -np.inf] + legend_handles = [] + xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$'] + ax.set_xlim(0, len(xticks)+1) + ax.set_xticklabels([''] + xticks + ['']) + ax.set_xlabel(r'BSM operator angle') + ylabel = r'${\rm log}_{10} \Lambda' + get_units(dim) + r'$' + ax.set_ylabel(ylabel) + for i_th in xrange(len(xticks)): + for i_frs, frs in enumerate(fix_sfr_mfr): + scale, evidences = angles_lim_array[i_dim][i_frs][i_th].T + null = evidences[0] + # TODO(shivesh): negative or not? + reduced_ev = -(evidences - null) + al = scale[reduced_ev > significance] + if len(al) > 0: + label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5]) + lim = al[0] + print 'frs, dim, lim = ', frs, dim, lim + if lim < yranges[0]: yranges[0] = lim + if lim > yranges[1]: yranges[1] = lim+4 + line = plt.Line2D( + (i_th+1-0.1, i_th+1+0.1), (lim, lim), lw=3, color=colour[i_frs], label=label + ) + ax.add_line(line) + if i_th == 0: legend_handles.append(line) + x_offset = i_frs*0.05 - 0.05 + ax.annotate( + s='', xy=(i_th+1+x_offset, lim), xytext=(i_th+1+x_offset, lim+3), + arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[i_frs]} + ) + 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', + title='dimension {0}'.format(dim)) + + 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) + + for of in outformat: + fig.savefig('./images/angles_limit_DIM{0}'.format(dim)+'.'+of, bbox_inches='tight', dpi=150) diff --git a/submitter/clean.sh b/submitter/clean.sh index edf0eea..a89f017 100755 --- a/submitter/clean.sh +++ b/submitter/clean.sh @@ -1,5 +1,7 @@ #!/bin/bash rm -f dagman_FR.submit.* +rm -f dagman_FR_*.submit.* rm -f logs/* rm -f metaouts/* +rm -f mnrun/* diff --git a/submitter/make_dag.py b/submitter/make_dag.py index 53878a2..78b7bff 100644 --- a/submitter/make_dag.py +++ b/submitter/make_dag.py @@ -19,10 +19,10 @@ fix_sfr_mfr = [ (1, 1, 1, 1, 2, 0), # (1, 1, 0, 1, 2, 0), # (1, 2, 0, 1, 2, 0), - # (1, 1, 1, 1, 0, 0), + (1, 1, 1, 1, 0, 0), # (1, 1, 0, 1, 0, 0), # (1, 0, 0, 1, 0, 0), - # (1, 1, 1, 0, 1, 0), + (1, 1, 1, 0, 1, 0), # (1, 1, 0, 0, 1, 0), # (0, 1, 0, 0, 1, 0), # (1, 2, 0, 0, 1, 0), @@ -31,17 +31,16 @@ fix_sfr_mfr = [ # MCMC run_mcmc = 'False' -burnin = 500 -nsteps = 2000 +burnin = 2500 +nsteps = 10000 nwalkers = 60 seed = 24 -threads = 1 +threads = 4 mcmc_seed_type = 'uniform' # FR dimension = [6] energy = [1e6] -likelihood = 'golemfit' no_bsm = 'False' sigma_ratio = ['0.01'] scale = "1E-20 1E-30" @@ -67,23 +66,30 @@ ast = 'p2_0' data = 'real' # Bayes Factor -run_bayes_factor = 'False' -run_angles_limit = 'True' -bayes_bins = 10 -bayes_live_points = 200 -bayes_tolerance = 0.01 -bayes_eval_bin = True # set to 'all' to run normally +run_bayes_factor = 'False' +run_angles_limit = 'False' +run_angles_correlation = 'True' +bayes_bins = 20 +bayes_live_points = 1000 +bayes_tolerance = 0.01 +bayes_eval_bin = 'None' # set to 'all' to run normally # Plot -plot_angles = 'False' -plot_elements = 'False' -plot_bayes = 'False' +plot_angles = 'False' +plot_elements = 'False' +plot_bayes = 'False' +plot_angles_limit = 'False' -outfile = 'dagman_FR_angles_limit.submit' +# outfile = 'dagman_FR.submit'.format(dimension[0]) +outfile = 'dagman_FR_angles_correlation_DIM{0}.submit'.format(dimension[0]) golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit' condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/submit.sub' -if bayes_eval_bin != 'all': b_runs = bayes_bins +if bayes_eval_bin != 'all': + if run_angles_correlation == 'True': + b_runs = bayes_bins**2 + else: + b_runs = bayes_bins else: b_runs = 1 with open(outfile, 'w') as f: @@ -99,6 +105,7 @@ with open(outfile, 'w') as f: outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(likelihood, dim, spectral_index) bayes_output = 'None' + angles_lim_output = 'None' for sig in sigma_ratio: print 'sigma', sig for frs in fix_sfr_mfr: @@ -108,6 +115,8 @@ with open(outfile, 'w') as f: bayes_output = outchains + '/bayes_factor/' if run_angles_limit == 'True': angles_lim_output = outchains + '/angles_limit/' + if run_angles_correlation == 'True': + angles_corr_output = outchains + '/angles_corr/' outchains += 'mcmc_chain' for r in range(b_runs): print 'run', r @@ -144,7 +153,6 @@ with open(outfile, 'w') as f: f.write('VARS\tjob{0}\tplot_elements="{1}"\n'.format(job_number, plot_elements)) f.write('VARS\tjob{0}\tseed="{1}"\n'.format(job_number, seed)) f.write('VARS\tjob{0}\tthreads="{1}"\n'.format(job_number, threads)) - f.write('VARS\tjob{0}\tlikelihood="{1}"\n'.format(job_number, likelihood)) f.write('VARS\tjob{0}\tmcmc_seed_type="{1}"\n'.format(job_number, mcmc_seed_type)) f.write('VARS\tjob{0}\tenergy_dependance="{1}"\n'.format(job_number, energy_dependance)) f.write('VARS\tjob{0}\tspectral_index="{1}"\n'.format(job_number, spectral_index)) @@ -161,6 +169,9 @@ with open(outfile, 'w') as f: f.write('VARS\tjob{0}\tbayes_eval_bin="{1}"\n'.format(job_number, r)) f.write('VARS\tjob{0}\trun_angles_limit="{1}"\n'.format(job_number, run_angles_limit)) f.write('VARS\tjob{0}\tangles_lim_output="{1}"\n'.format(job_number, angles_lim_output)) + f.write('VARS\tjob{0}\tplot_angles_limit="{1}"\n'.format(job_number, plot_angles_limit)) + f.write('VARS\tjob{0}\trun_angles_correlation="{1}"\n'.format(job_number, run_angles_correlation)) + f.write('VARS\tjob{0}\tangles_corr_output="{1}"\n'.format(job_number, angles_corr_output)) job_number += 1 for frs in full_scan_mfr: @@ -170,6 +181,8 @@ with open(outfile, 'w') as f: bayes_output = outchains + '/bayes_factor/' if run_angles_limit == 'True': angles_lim_output = outchains + '/angles_limit/' + if run_angles_correlation == 'True': + angles_corr_output = outchains + '/angles_corr/' outchains += 'mcmc_chain' for r in range(b_runs): print 'run', r @@ -206,7 +219,6 @@ with open(outfile, 'w') as f: f.write('VARS\tjob{0}\tplot_elements="{1}"\n'.format(job_number, plot_elements)) f.write('VARS\tjob{0}\tseed="{1}"\n'.format(job_number, seed)) f.write('VARS\tjob{0}\tthreads="{1}"\n'.format(job_number, threads)) - f.write('VARS\tjob{0}\tlikelihood="{1}"\n'.format(job_number, likelihood)) f.write('VARS\tjob{0}\tmcmc_seed_type="{1}"\n'.format(job_number, mcmc_seed_type)) f.write('VARS\tjob{0}\tenergy_dependance="{1}"\n'.format(job_number, energy_dependance)) f.write('VARS\tjob{0}\tspectral_index="{1}"\n'.format(job_number, spectral_index)) @@ -223,4 +235,7 @@ with open(outfile, 'w') as f: f.write('VARS\tjob{0}\tbayes_eval_bin="{1}"\n'.format(job_number, r)) f.write('VARS\tjob{0}\trun_angles_limit="{1}"\n'.format(job_number, run_angles_limit)) f.write('VARS\tjob{0}\tangles_lim_output="{1}"\n'.format(job_number, angles_lim_output)) + f.write('VARS\tjob{0}\tplot_angles_limit="{1}"\n'.format(job_number, plot_angles_limit)) + f.write('VARS\tjob{0}\trun_angles_correlation="{1}"\n'.format(job_number, run_angles_correlation)) + f.write('VARS\tjob{0}\tangles_corr_output="{1}"\n'.format(job_number, angles_corr_output)) job_number += 1 diff --git a/submitter/submit.sub b/submitter/submit.sub index d232097..810ffa0 100644 --- a/submitter/submit.sub +++ b/submitter/submit.sub @@ -1,5 +1,5 @@ Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/fr.py -Arguments = "--ast $(ast) --astroDeltaGamma $(astroDeltaGamma) --astroNorm $(astroNorm) --burnin $(burnin) --convNorm $(convNorm) --data $(data) --dimension $(dimension) --energy $(energy) --fix-mixing $(fix_mixing) --fix-scale $(fix_scale) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --muonNorm $(muonNorm) --no-bsm $(no_bsm) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --promptNorm $(promptNorm) --run-mcmc $(run_mcmc) --scale $(scale) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --likelihood $(likelihood) --mcmc-seed-type $(mcmc_seed_type) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning_0) $(binning_1) $(binning_2) --fix-mixing-almost $(fix_mixing_almost) --run-bayes-factor $(run_bayes_factor) --bayes-bins $(bayes_bins) --bayes-output $(bayes_output) --bayes-live-points $(bayes_live_points) --plot-bayes $(plot_bayes) --bayes-tolerance $(bayes_tolerance) --bayes-eval-bin $(bayes_eval_bin) --run-angles-limit $(run_angles_limit) --angles-lim-out $(angles_lim_output)" +Arguments = "--ast $(ast) --astroDeltaGamma $(astroDeltaGamma) --astroNorm $(astroNorm) --burnin $(burnin) --convNorm $(convNorm) --data $(data) --dimension $(dimension) --energy $(energy) --fix-mixing $(fix_mixing) --fix-scale $(fix_scale) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --muonNorm $(muonNorm) --no-bsm $(no_bsm) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --promptNorm $(promptNorm) --run-mcmc $(run_mcmc) --scale $(scale) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --mcmc-seed-type $(mcmc_seed_type) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning_0) $(binning_1) $(binning_2) --fix-mixing-almost $(fix_mixing_almost) --run-bayes-factor $(run_bayes_factor) --bayes-bins $(bayes_bins) --bayes-output $(bayes_output) --bayes-live-points $(bayes_live_points) --plot-bayes $(plot_bayes) --bayes-tolerance $(bayes_tolerance) --bayes-eval-bin $(bayes_eval_bin) --run-angles-limit $(run_angles_limit) --angles-lim-out $(angles_lim_output) --plot-angles-limit $(plot_angles_limit) --run-angles-correlation $(run_angles_correlation) --angles-corr-output $(angles_corr_output)" # All logs will go to a single file log = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log @@ -14,9 +14,9 @@ getenv = True # but do not try to copy outputs back (see: https://htcondor-wiki.cs.wisc.edu/index.cgi/tktview?tn=3081) # +TransferOutput="" -transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/ +Transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/ -request_memory = 10GB +request_memory = 1GB 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/utils/plot.py b/utils/plot.py index 0d02c2a..0b82675 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -264,7 +264,7 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): g.export(outfile+'_elements.'+of) -def bayes_factor_plot(dirname, outfile, outformat, args, xlim): +def bayes_factor_plot(dirname, outfile, outformat, args): """Make Bayes factor plot.""" if not args.plot_bayes: return print "Making Bayes Factor plot" @@ -287,7 +287,7 @@ def bayes_factor_plot(dirname, outfile, outformat, args, xlim): fig = plt.figure(figsize=(7, 5)) ax = fig.add_subplot(111) - ax.set_xlim(xlim) + ax.set_xlim(np.log10(args.scale_region)) ax.set_xlabel(r'${\rm log}_{10} \Lambda ' + get_units(args.dimension) +r'$') ax.set_ylabel(r'Bayes Factor') |