aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xsens.py164
-rw-r--r--submitter/make_dag.py241
-rw-r--r--submitter/mcmc_dag.py116
-rw-r--r--submitter/mcmc_submit.sub39
-rw-r--r--submitter/sens_dag.py134
-rw-r--r--submitter/sens_submit.sub39
-rw-r--r--submitter/submit.sub41
-rw-r--r--utils/enums.py13
-rw-r--r--utils/fr.py49
-rw-r--r--utils/gf.py22
-rw-r--r--utils/likelihood.py30
-rw-r--r--utils/misc.py5
-rw-r--r--utils/multinest.py37
-rw-r--r--utils/param.py9
-rw-r--r--utils/plot.py9
15 files changed, 515 insertions, 433 deletions
diff --git a/sens.py b/sens.py
index 980ef3f..50fe0c8 100755
--- a/sens.py
+++ b/sens.py
@@ -25,6 +25,7 @@ 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.misc import DTYPE
from utils.param import Param, ParamSet, get_paramsets
from utils import multinest as mn_utils
@@ -34,11 +35,13 @@ 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.29, 0.31], 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.998, 1.0], 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.46, 0.61], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag),
- Param(name='dcp', value=4.08404, seed=[0, 2*np.pi], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag),
+ 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
@@ -76,9 +79,9 @@ def process_args(args):
raise NotImplementedError(
'--fix-mixing and --fix-mixing-almost cannot be used together'
)
- if args.mn_run and args.fix_scale:
+ if args.sens_run and args.fix_scale:
raise NotImplementedError(
- '--mn-run and --fix-scale cannot be used together'
+ '--sens-run and --fix-scale cannot be used together'
)
args.measured_ratio = fr_utils.normalise_fr(args.measured_ratio)
@@ -94,13 +97,13 @@ def process_args(args):
args.scale = fr_utils.estimate_scale(args)
args.scale_region = (args.scale/args.scale_region, args.scale*args.scale_region)
- if args.mn_eval_bin.lower() == 'all':
- args.mn_eval_bin = None
+ if args.sens_eval_bin.lower() == 'all':
+ args.sens_eval_bin = None
else:
- args.mn_eval_bin = int(args.mn_eval_bin)
+ args.sens_eval_bin = int(args.sens_eval_bin)
if args.stat_method is StatCateg.FREQUENTIST and \
- args.likelihood is Likelihood.GOLEMFIT::
+ args.likelihood is Likelihood.GOLEMFIT:
args.likelihood = Likelihood.GF_FREQ
@@ -123,6 +126,10 @@ def parse_args(args=None):
help='Path to output chains'
)
parser.add_argument(
+ '--sens-run', type=misc_utils.parse_bool, default='True',
+ help='Generate sensitivities'
+ )
+ parser.add_argument(
'--run-method', default='full',
type=partial(misc_utils.enum_parse, c=SensitivityCateg),
choices=SensitivityCateg,
@@ -134,7 +141,15 @@ def parse_args(args=None):
help='Statistical method to employ'
)
parser.add_argument(
- '--plot-statistic', type=parse_bool, default='False',
+ '--sens-bins', type=int, default=10,
+ help='Number of bins for the Bayes factor plot'
+ )
+ parser.add_argument(
+ '--sens-eval-bin', type=str, default='all',
+ help='Which bin to evalaute for Bayes factor plot'
+ )
+ parser.add_argument(
+ '--plot-statistic', type=misc_utils.parse_bool, default='False',
help='Plot MultiNest evidence or LLH value'
)
fr_utils.fr_argparse(parser)
@@ -161,8 +176,10 @@ def main():
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 is SensitivityCateg.FIXED_ANGLE or \
- args.run_method is SensitivityCateg.CORR_ANGLE:
+ 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):
@@ -171,55 +188,66 @@ def main():
if mmangles[x].name != prms.name])
)
+ scan_scales = np.linspace(
+ 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
+
+ if args.run_method is SensitivityCateg.CORR_ANGLE:
+ scan_angles = np.linspace(0+1e-9, 1-1e-9, eval_dim)
+ else: scan_angles = np.array([0])
+ print 'scan_scales', scan_scales
+ print 'scan_angles', scan_angles
+
out = args.outfile+'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
+ misc_utils.gen_identifier(args)
- if args.mn_run:
- if args.likelihood is Likelihood.GOLEMFIT:
+ if args.sens_run:
+ if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
fitter = gf_utils.setup_fitter(args, asimov_paramset)
- if args.StatCateg is StatCateg.FREQUENTIST:
- fitter.SetFitParametersFlag(gf_utils.fit_flags(llh_paramset)
+ if args.stat_method is StatCateg.FREQUENTIST:
+ flags, gf_nuisance = gf_utils.fit_flags(llh_paramset)
+ llh_paramset = llh_paramset.remove_params(gf_nuisance)
+ asimov_paramset = asimov_paramset.remove_params(gf_nuisance)
+ st_paramset_arr = [x.remove_params(gf_nuisance)
+ for x in st_paramset_arr]
+ fitter.SetFitParametersFlag(flags)
else: fitter = None
- scan_scales = np.linspace(
- np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.mn_bins
- )
- if args.run_method is SensitivityCateg.CORR_ANGLE:
- scan_angles = np.linspace(0, 1, eval_dim)
- else: scan_angles = np.array([0])
- print 'scan_scales', scan_scales
- print 'scan_angles', scan_angles
-
- if args.mn_eval_bin is None:
- eval_dim = args.mn_bins
- else: eval_dim = 1
-
if args.run_method is SensitivityCateg.FULL:
statistic_arr = np.full((eval_dim, 2), np.nan)
elif args.run_method is SensitivityCateg.FIXED_ANGLE:
statistic_arr = np.full((len(st_paramset_arr), eval_dim, 2), np.nan)
elif args.run_method is SensitivityCateg.CORR_ANGLE:
- statistic_arr = np.full((len(st_paramset_arr), eval_dim, eval_dim, 3), np.nan)
-
+ statistic_arr = np.full(
+ (len(st_paramset_arr), eval_dim, eval_dim, 3), np.nan
+ )
- for idx_scen, mn_paramset in enumerate(st_paramset_arr):
+ for idx_scen, sens_paramset in enumerate(st_paramset_arr):
print '|||| SCENARIO = {0}'.format(idx_scen)
- for x in mmangles: x.value = 0.
- if args.run_method is SensitivityCateg.FIXED_ANGLE:
+ 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 idx_scen == 0 or idx_scen == 2:
- mmangles[idx_scen].value = np.sin(np.pi/2.)**2
+ mmangles[idx_scen].value = np.sin(np.pi/4., dtype=DTYPE)**2
"""s_12^2 or s_23^2"""
elif idx_scen == 1:
- mmangles[idx_scen].value = np.cos(np.pi/2.)**4
+ mmangles[idx_scen].value = np.cos(np.pi/4., dtype=DTYPE)**4
"""c_13^4"""
for idx_an, an in enumerate(scan_angles):
- if args.run_method is SensitivityCateg.CORR_ANGLE:
+ if args.run_method in [SensitivityCateg.CORR_ANGLE,
+ SensitivityCateg.CORR_ONE_ANGLE]:
print '|||| ANGLE = {0:<04.2}'.format(an)
- mmangles[idx_an].value = an
+ if SensitivityCateg.CORR_ANGLE:
+ for x in mmangles: x.value = 0.+1e-9
+ mmangles[idx_scen].value = an
for idx_sc, sc in enumerate(scan_scales):
- if args.mn_eval_bin is not None:
- if idx_sc == args.mn_eval_bin:
+ 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))
@@ -231,7 +259,7 @@ def main():
if args.stat_method is StatCateg.BAYESIAN:
try:
stat = mn_utils.mn_evidence(
- mn_paramset = mn_paramset,
+ mn_paramset = sens_paramset,
llh_paramset = llh_paramset,
asimov_paramset = asimov_paramset,
args = args,
@@ -242,26 +270,30 @@ def main():
continue
print '## Evidence = {0}'.format(stat)
elif args.stat_method is StatCateg.FREQUENTIST:
- llh_paramset_freq = [x for parm in llh_paramset if
- x.name not in asimov_paramset.names]
def fn(x):
- for idx, parm in enumerate(llh_paramset_freq):
- parm.value = x[idx]
- theta = llh_paramset_freq.values
- try:
- llh = llh_utils.ln_prob(
- theta=theta, args=args, asimov_paramset=asimov_paramset,
- mcmc_paramset=mcmc_paramset_freq, fitter=fitter
- )
- except:
- print 'Failed run, continuing'
- return np.inf
+ pranges = sens_paramset.seeds
+ for i, name in enumerate(sens_paramset.names):
+ llh_paramset[name].value = \
+ (pranges[i][1]-pranges[i][0])*x[i] + pranges[i][0]
+ theta = llh_paramset.values
+ print 'llh_paramset', llh_paramset
+ llh = llh_utils.ln_prob(
+ theta=theta, args=args, asimov_paramset=asimov_paramset,
+ llh_paramset=llh_paramset, fitter=fitter
+ )
+ print 'llh', llh
return -llh
-
- x0 = np.average(llh_paramset_freq.values, axis=1)
- res = minimize(fun=fn, x0=x0, method='L-BFGS-B',
- bounds=llh_paramset.seed, fitter=fitter)
+
+ n_params = len(sens_paramset)
+ x0 = np.full(n_params, 0.7)
+ bounds = np.dstack([np.zeros(n_params), np.ones(n_params)])[0]
+ try:
+ res = minimize(fun=fn, x0=x0, method='L-BFGS-B', bounds=bounds)
+ except AssertionError:
+ print 'Failed run, continuing'
+ continue
stat = -fn(res.x)
+ 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:
@@ -271,10 +303,10 @@ def main():
misc_utils.make_dir(out)
print 'Saving to {0}'.format(out+'.npy')
- np.save(out+'.npy', np.array(statistic_arr))
+ np.save(out+'.npy', statistic_arr)
if args.plot_statistic:
- if args.mn_run: raw = statistic_arr
+ if args.sens_run: raw = statistic_arr
else: raw = np.load(out+'.npy')
data = ma.masked_invalid(raw, 0)
@@ -292,9 +324,9 @@ def main():
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{1}{2}$'
- elif idx_scen == 1: label = r'$\mathcal{O}_{13}=\frac{1}{2}$'
- elif idx_scen == 2: label = r'$\mathcal{O}_{23}=\frac{1}{2}$'
+ 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[idx_scen],
outfile = outfile,
@@ -313,14 +345,14 @@ def main():
for idx_an, an in enumerate(scan_angles):
print '|||| ANGLE = {0:<04.2}'.format(an)
outfile = basescenoutfile + '_ANGLE{0}'.format(idx_an)
- label += r'{0:<04.2}$'.format(an)
+ _label = label + r'{0:<04.2}$'.format(an)
plot_utils.plot_statistic(
data = data[idx_scen][idx_an][:,1:],
outfile = outfile,
outformat = ['png'],
args = args,
scale_param = scale,
- label = label
+ label = _label
)
diff --git a/submitter/make_dag.py b/submitter/make_dag.py
deleted file mode 100644
index 0e41c9a..0000000
--- a/submitter/make_dag.py
+++ /dev/null
@@ -1,241 +0,0 @@
-#! /usr/bin/env python
-
-import os
-import numpy as np
-
-a_fr = (1, 2, 0)
-b_fr = (1, 0, 0)
-c_fr = (0, 1, 0)
-d_fr = (0, 0, 1)
-e_fr = (1, 1, 1)
-f_fr = (2, 1, 0)
-g_fr = (1, 1, 0)
-
-full_scan_mfr = [
- # (1, 1, 1), (1, 1, 0)
-]
-
-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, 0, 1, 2, 0),
- # (1, 1, 0, 1, 0, 0),
- # (1, 1, 0, 0, 1, 0),
- # (1, 0, 0, 1, 0, 0),
- # (0, 1, 0, 0, 1, 0),
- # (1, 2, 0, 1, 2, 0),
- # (1, 2, 0, 0, 1, 0),
-]
-
-# MCMC
-run_mcmc = 'True'
-burnin = 2500
-nsteps = 10000
-nwalkers = 60
-seed = 'None'
-threads = 1
-mcmc_seed_type = 'uniform'
-
-# FR
-dimension = [3, 6]
-energy = [1e6]
-no_bsm = 'False'
-sigma_ratio = ['0.01']
-scale = "1E-20 1E-30"
-scale_region = "1E10"
-energy_dependance = 'spectral'
-spectral_index = -2
-binning = [1e4, 1e7, 5]
-fix_mixing = 'False'
-fix_mixing_almost = 'False'
-
-# Likelihood
-likelihood = 'gaussian'
-
-# Nuisance
-convNorm = 1.
-promptNorm = 0.
-muonNorm = 1.
-astroNorm = 6.9
-astroDeltaGamma = 2.5
-
-# GolemFit
-ast = 'p2_0'
-data = 'real'
-
-# Bayes Factor
-run_bayes_factor = 'False'
-run_angles_limit = 'False'
-run_angles_correlation = 'False'
-bayes_bins = 100
-bayes_live_points = 3000
-bayes_tolerance = 0.01
-bayes_eval_bin = 'all' # set to 'all' to run normally
-
-# Plot
-plot_angles = 'True'
-plot_elements = 'False'
-plot_bayes = 'False'
-plot_angles_limit = 'False'
-
-outfile = 'dagman_FR.submit'
-golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
-condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/submit.sub'
-
-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:
- job_number = 1
- for dim in dimension:
- print 'dimension', dim
- for en in energy:
- print 'energy {0:.0E}'.format(en)
-
- 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:
- print 'sigma', sig
- for frs in fix_sfr_mfr:
- print frs
- outchains = outchain_head + '/fix_ifr/{0}/'.format(str(sig).replace('.', '_'))
- if run_bayes_factor == 'True':
- 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
- f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
- f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
- f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
- f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
- f.write('VARS\tjob{0}\tsigma_ratio="{1}"\n'.format(job_number, sig))
- f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, 'True'))
- 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]))
- f.write('VARS\tjob{0}\tfix_scale="{1}"\n'.format(job_number, 'False'))
- f.write('VARS\tjob{0}\tscale="{1}"\n'.format(job_number, 0))
- f.write('VARS\tjob{0}\tscale_region="{1}"\n'.format(job_number, scale_region))
- f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
- f.write('VARS\tjob{0}\tenergy="{1}"\n'.format(job_number, en))
- f.write('VARS\tjob{0}\tlikelihood="{1}"\n'.format(job_number, likelihood))
- f.write('VARS\tjob{0}\tburnin="{1}"\n'.format(job_number, burnin))
- f.write('VARS\tjob{0}\tnwalkers="{1}"\n'.format(job_number, nwalkers))
- f.write('VARS\tjob{0}\tnsteps="{1}"\n'.format(job_number, nsteps))
- f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
- f.write('VARS\tjob{0}\tfix_mixing="{1}"\n'.format(job_number, fix_mixing))
- f.write('VARS\tjob{0}\tno_bsm="{1}"\n'.format(job_number, no_bsm))
- f.write('VARS\tjob{0}\trun_mcmc="{1}"\n'.format(job_number, run_mcmc))
- f.write('VARS\tjob{0}\tastroDeltaGamma="{1}"\n'.format(job_number, astroDeltaGamma))
- f.write('VARS\tjob{0}\tastroNorm="{1}"\n'.format(job_number, astroNorm))
- f.write('VARS\tjob{0}\tconvNorm="{1}"\n'.format(job_number, convNorm))
- f.write('VARS\tjob{0}\tmuonNorm="{1}"\n'.format(job_number, muonNorm))
- f.write('VARS\tjob{0}\tpromptNorm="{1}"\n'.format(job_number, promptNorm))
- f.write('VARS\tjob{0}\tdata="{1}"\n'.format(job_number, data))
- f.write('VARS\tjob{0}\tast="{1}"\n'.format(job_number, ast))
- f.write('VARS\tjob{0}\tplot_angles="{1}"\n'.format(job_number, plot_angles))
- 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}\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))
- f.write('VARS\tjob{0}\tbinning_0="{1}"\n'.format(job_number, binning[0]))
- f.write('VARS\tjob{0}\tbinning_1="{1}"\n'.format(job_number, binning[1]))
- f.write('VARS\tjob{0}\tbinning_2="{1}"\n'.format(job_number, binning[2]))
- f.write('VARS\tjob{0}\tfix_mixing_almost="{1}"\n'.format(job_number, fix_mixing_almost))
- # f.write('VARS\tjob{0}\trun_bayes_factor="{1}"\n'.format(job_number, run_bayes_factor))
- # f.write('VARS\tjob{0}\tbayes_bins="{1}"\n'.format(job_number, bayes_bins))
- # f.write('VARS\tjob{0}\tbayes_output="{1}"\n'.format(job_number, bayes_output))
- # f.write('VARS\tjob{0}\tbayes_live_points="{1}"\n'.format(job_number, bayes_live_points))
- # f.write('VARS\tjob{0}\tbayes_tolerance="{1}"\n'.format(job_number, bayes_tolerance))
- # f.write('VARS\tjob{0}\tplot_bayes="{1}"\n'.format(job_number, plot_bayes))
- # 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:
- print frs
- outchains = outchain_head + '/full_scan/{0}'.format(str(sig).replace('.', '_'))
- if run_bayes_factor == 'True':
- 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
- f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
- f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
- f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
- f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
- f.write('VARS\tjob{0}\tsigma_ratio="{1}"\n'.format(job_number, sig))
- f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, 'False'))
- 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))
- f.write('VARS\tjob{0}\tfix_scale="{1}"\n'.format(job_number, 'False'))
- f.write('VARS\tjob{0}\tscale="{1}"\n'.format(job_number, 0))
- f.write('VARS\tjob{0}\tscale_region="{1}"\n'.format(job_number, scale_region))
- f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
- f.write('VARS\tjob{0}\tenergy="{1}"\n'.format(job_number, en))
- f.write('VARS\tjob{0}\tlikelihood="{1}"\n'.format(job_number, likelihood))
- f.write('VARS\tjob{0}\tburnin="{1}"\n'.format(job_number, burnin))
- f.write('VARS\tjob{0}\tnwalkers="{1}"\n'.format(job_number, nwalkers))
- f.write('VARS\tjob{0}\tnsteps="{1}"\n'.format(job_number, nsteps))
- f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
- f.write('VARS\tjob{0}\tfix_mixing="{1}"\n'.format(job_number, fix_mixing))
- f.write('VARS\tjob{0}\tno_bsm="{1}"\n'.format(job_number, no_bsm))
- f.write('VARS\tjob{0}\trun_mcmc="{1}"\n'.format(job_number, run_mcmc))
- f.write('VARS\tjob{0}\tastroDeltaGamma="{1}"\n'.format(job_number, astroDeltaGamma))
- f.write('VARS\tjob{0}\tastroNorm="{1}"\n'.format(job_number, astroNorm))
- f.write('VARS\tjob{0}\tconvNorm="{1}"\n'.format(job_number, convNorm))
- f.write('VARS\tjob{0}\tmuonNorm="{1}"\n'.format(job_number, muonNorm))
- f.write('VARS\tjob{0}\tpromptNorm="{1}"\n'.format(job_number, promptNorm))
- f.write('VARS\tjob{0}\tdata="{1}"\n'.format(job_number, data))
- f.write('VARS\tjob{0}\tast="{1}"\n'.format(job_number, ast))
- f.write('VARS\tjob{0}\tplot_angles="{1}"\n'.format(job_number, plot_angles))
- 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}\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))
- f.write('VARS\tjob{0}\tbinning_0="{1}"\n'.format(job_number, binning[0]))
- f.write('VARS\tjob{0}\tbinning_1="{1}"\n'.format(job_number, binning[1]))
- f.write('VARS\tjob{0}\tbinning_2="{1}"\n'.format(job_number, binning[2]))
- f.write('VARS\tjob{0}\tfix_mixing_almost="{1}"\n'.format(job_number, fix_mixing_almost))
- f.write('VARS\tjob{0}\trun_bayes_factor="{1}"\n'.format(job_number, run_bayes_factor))
- f.write('VARS\tjob{0}\tbayes_bins="{1}"\n'.format(job_number, bayes_bins))
- f.write('VARS\tjob{0}\tbayes_output="{1}"\n'.format(job_number, bayes_output))
- f.write('VARS\tjob{0}\tbayes_live_points="{1}"\n'.format(job_number, bayes_live_points))
- f.write('VARS\tjob{0}\tbayes_tolerance="{1}"\n'.format(job_number, bayes_tolerance))
- f.write('VARS\tjob{0}\tplot_bayes="{1}"\n'.format(job_number, plot_bayes))
- 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/mcmc_dag.py b/submitter/mcmc_dag.py
new file mode 100644
index 0000000..79f9b6d
--- /dev/null
+++ b/submitter/mcmc_dag.py
@@ -0,0 +1,116 @@
+#! /usr/bin/env python
+
+import os
+import numpy as np
+
+full_scan_mfr = [
+ # (1, 1, 1), (1, 0, 0)
+]
+
+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, 0, 1, 2, 0),
+ # (1, 1, 0, 1, 0, 0),
+ # (1, 1, 0, 0, 1, 0),
+ # (1, 0, 0, 1, 0, 0),
+ # (0, 1, 0, 0, 1, 0),
+ # (1, 2, 0, 1, 2, 0),
+ # (1, 2, 0, 0, 1, 0),
+]
+
+GLOBAL_PARAMS = {}
+
+# MCMC
+GLOBAL_PARAMS.update(dict(
+ run_mcmc = 'True',
+ burnin = 2500,
+ nsteps = 10000,
+ nwalkers = 60,
+ seed = 'None',
+ mcmc_seed_type = 'uniform'
+))
+
+# FR
+dimension = [3, 6]
+GLOBAL_PARAMS.update(dict(
+ threads = 1,
+ binning = '1e4 1e7 5',
+ no_bsm = 'False',
+ scale_region = "1E10",
+ energy_dependance = 'spectral',
+ spectral_index = -2,
+ fix_mixing = 'False',
+ fix_mixing_almost = 'False'
+))
+
+# Likelihood
+GLOBAL_PARAMS.update(dict(
+ likelihood = 'gaussian',
+ sigma_ratio = '0.01'
+))
+
+# GolemFit
+GLOBAL_PARAMS.update(dict(
+ ast = 'p2_0',
+ data = 'real'
+))
+
+# Plot
+GLOBAL_PARAMS.update(dict(
+ plot_angles = 'True',
+ plot_elements = 'False',
+))
+
+outfile = 'dagman_FR_MCMC.submit'
+golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
+condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/mcmc_submit.sub'
+
+with open(outfile, 'w') as f:
+ job_number = 1
+ for dim in dimension:
+ print 'dimension', dim
+ outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(
+ GLOBAL_PARAMS['likelihood'], dim, GLOBAL_PARAMS['spectral_index']
+ )
+ for frs in fix_sfr_mfr:
+ print 'frs', frs
+ outchains = outchain_head + '/fix_ifr/'
+ if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
+ outchains += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+ outchains += 'mcmc_chain'
+ f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+ f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+ f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+ f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+ f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+ f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, True))
+ 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))
+ job_number += 1
+
+ for frs in full_scan_mfr:
+ print 'frs', frs
+ outchains = outchain_head + '/full/'
+ if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
+ outchains += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+ outchains += 'mcmc_chain'
+ f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+ f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+ f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+ f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+ f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+ f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, False))
+ 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))
+ job_number += 1
diff --git a/submitter/mcmc_submit.sub b/submitter/mcmc_submit.sub
new file mode 100644
index 0000000..2032cb6
--- /dev/null
+++ b/submitter/mcmc_submit.sub
@@ -0,0 +1,39 @@
+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)"
+
+# All logs will go to a single file
+log = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
+output = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).out
+error = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).err
+
+getenv = True
+# environment = "X509_USER_PROXY=x509up_u14830"
+
+# Stage user cert to the node (Gridftp-Users is already on CVMFS)
+# transfer_input_files = /tmp/x509up_u14830
+
+# 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/
+
+request_memory = 1GB
+request_cpus = 1
+
+Universe = vanilla
+Notification = never
+
+# run on both SL5 and 6
+# +WantRHEL6 = True
+# +WantSLC6 = False
+
+# # run on OSG
+# +WantGlidein = True
+
+# +TransferOutput=""
+
+# Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
+# Requirements = IS_GLIDEIN
+# Requirements = (OpSysMajorVer =?= 6)
+
+# GO!
+queue
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
new file mode 100644
index 0000000..2cecfe9
--- /dev/null
+++ b/submitter/sens_dag.py
@@ -0,0 +1,134 @@
+#! /usr/bin/env python
+
+import os
+import numpy as np
+
+full_scan_mfr = [
+ # (1, 1, 1), (1, 2, 0)
+]
+
+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, 0, 1, 2, 0),
+ # (1, 1, 0, 1, 0, 0),
+ # (1, 1, 0, 0, 1, 0),
+ # (1, 0, 0, 1, 0, 0),
+ # (0, 1, 0, 0, 1, 0),
+ # (1, 2, 0, 1, 2, 0),
+ # (1, 2, 0, 0, 1, 0),
+]
+
+GLOBAL_PARAMS = {}
+
+# Bayes Factor
+GLOBAL_PARAMS.update(dict(
+ sens_run = 'True',
+ run_method = 'full',
+ stat_method = 'bayesian',
+ sens_bins = 10,
+ sens_eval_bin = 'all' # set to 'all' to run normally
+))
+
+# MultiNest
+GLOBAL_PARAMS.update(dict(
+ mn_live_points = 400,
+ mn_tolerance = 0.01,
+ mn_output = './mnrun'
+))
+
+# FR
+dimension = [3, 6]
+GLOBAL_PARAMS.update(dict(
+ threads = 1,
+ binning = '1e4 1e7 5',
+ no_bsm = 'False',
+ scale_region = "1E10",
+ energy_dependance = 'spectral',
+ spectral_index = -2,
+ fix_mixing = 'False',
+ fix_mixing_almost = 'False'
+))
+
+# Likelihood
+GLOBAL_PARAMS.update(dict(
+ likelihood = 'gaussian',
+ sigma_ratio = '0.01'
+))
+
+# GolemFit
+GLOBAL_PARAMS.update(dict(
+ ast = 'p2_0',
+ data = 'real'
+))
+
+# Plot
+GLOBAL_PARAMS.update(dict(
+ plot_statistic = 'True'
+))
+
+outfile = 'dagman_FR_SENS.submit'
+golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
+condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
+
+if GLOBAL_PARAMS['sens_eval_bin'].lower() != 'all':
+ if GLOBAL_PARAMS['run_method'].lower() == 'corr_angle':
+ sens_runs = GLOBAL_PARAMS['sens_bins']**2
+ else:
+ sens_runs = GLOBAL_PARAMS['sens_bins']
+else: sens_runs = 1
+
+with open(outfile, 'w') as f:
+ job_number = 1
+ for dim in dimension:
+ print 'dimension', dim
+ outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(
+ GLOBAL_PARAMS['likelihood'], dim, GLOBAL_PARAMS['spectral_index']
+ )
+ for frs in fix_sfr_mfr:
+ print 'frs', frs
+ output = outchain_head + '/fix_ifr/'
+ if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
+ output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+ output += 'fr_stat'
+ for r in xrange(sens_runs):
+ print 'run', r
+ f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+ f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+ f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+ f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+ f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+ f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, True))
+ 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]))
+ 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))
+ job_number += 1
+
+ for frs in full_scan_mfr:
+ print 'frs', frs
+ output = outchain_head + '/full/'
+ if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
+ output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+ output += 'fr_stat'
+ for r in xrange(sens_runs):
+ print 'run', r
+ f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+ f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+ f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+ f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+ f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+ f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, False))
+ 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))
+ 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))
+ job_number += 1
diff --git a/submitter/sens_submit.sub b/submitter/sens_submit.sub
new file mode 100644
index 0000000..1a02608
--- /dev/null
+++ b/submitter/sens_submit.sub
@@ -0,0 +1,39 @@
+Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/fr.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
+log = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
+output = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).out
+error = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).err
+
+getenv = True
+# environment = "X509_USER_PROXY=x509up_u14830"
+
+# Stage user cert to the node (Gridftp-Users is already on CVMFS)
+# transfer_input_files = /tmp/x509up_u14830
+
+# 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/
+
+request_memory = 1GB
+request_cpus = 1
+
+Universe = vanilla
+Notification = never
+
+# run on both SL5 and 6
+# +WantRHEL6 = True
+# +WantSLC6 = False
+
+# # run on OSG
+# +WantGlidein = True
+
+# +TransferOutput=""
+
+# Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
+# Requirements = IS_GLIDEIN
+# Requirements = (OpSysMajorVer =?= 6)
+
+# GO!
+queue
diff --git a/submitter/submit.sub b/submitter/submit.sub
deleted file mode 100644
index 8e3c5af..0000000
--- a/submitter/submit.sub
+++ /dev/null
@@ -1,41 +0,0 @@
-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) --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
-output = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).out
-error = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).err
-
-getenv = True
-# environment = "X509_USER_PROXY=x509up_u14830"
-
-# Stage user cert to the node (Gridftp-Users is already on CVMFS)
-# transfer_input_files = /tmp/x509up_u14830
-
-# 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/
-
-request_memory = 1GB
-request_cpus = 4
-
-Universe = vanilla
-Notification = never
-
-# run on both SL5 and 6
-# +WantRHEL6 = True
-# +WantSLC6 = False
-
-# # run on OSG
-# +WantGlidein = True
-
-# +TransferOutput=""
-
-# Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
-# Requirements = IS_GLIDEIN
-# Requirements = (OpSysMajorVer =?= 6)
-
-# GO!
-queue
diff --git a/utils/enums.py b/utils/enums.py
index 8a9e868..b80b712 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -39,8 +39,9 @@ class ParamTag(Enum):
class PriorsCateg(Enum):
- UNIFORM = 1
- GAUSSIAN = 2
+ UNIFORM = 1
+ GAUSSIAN = 2
+ HALFGAUSS = 3
class MCMCSeedType(Enum):
@@ -49,9 +50,11 @@ class MCMCSeedType(Enum):
class SensitivityCateg(Enum):
- FULL = 1
- FIXED_ANGLE = 2
- CORR_ANGLE = 3
+ FULL = 1
+ FIXED_ANGLE = 2
+ CORR_ANGLE = 3
+ FIXED_ONE_ANGLE = 4
+ CORR_ONE_ANGLE = 5
class StatCateg(Enum):
diff --git a/utils/fr.py b/utils/fr.py
index 342a848..a82e081 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -13,9 +13,9 @@ import argparse
from functools import partial
import numpy as np
-from scipy import linalg
from utils.enums import EnergyDependance
+from utils.misc import DTYPE, CDTYPE, PI
from utils.misc import enum_parse, parse_bool
@@ -44,11 +44,11 @@ def angles_to_fr(src_angles):
"""
sphi4, c2psi = src_angles
- psi = (0.5)*np.arccos(c2psi)
+ psi = (0.5)*np.arccos(c2psi, dtype=DTYPE)
- sphi2 = np.sqrt(sphi4)
+ sphi2 = np.sqrt(sphi4, dtype=DTYPE)
cphi2 = 1. - sphi2
- spsi2 = np.sin(psi)**2
+ spsi2 = np.sin(psi, dtype=DTYPE)**2
cspi2 = 1. - spsi2
x = sphi2*cspi2
@@ -80,16 +80,16 @@ def angles_to_u(bsm_angles):
"""
s12_2, c13_4, s23_2, dcp = bsm_angles
- dcp = np.complex128(dcp)
+ dcp = CDTYPE(dcp)
c12_2 = 1. - s12_2
- c13_2 = np.sqrt(c13_4)
+ c13_2 = np.sqrt(c13_4, dtype=DTYPE)
s13_2 = 1. - c13_2
c23_2 = 1. - s23_2
- t12 = np.arcsin(np.sqrt(s12_2))
- t13 = np.arccos(np.sqrt(c13_2))
- t23 = np.arcsin(np.sqrt(s23_2))
+ t12 = np.arcsin(np.sqrt(s12_2, dtype=DTYPE))
+ t13 = np.arccos(np.sqrt(c13_2, dtype=DTYPE))
+ t23 = np.arcsin(np.sqrt(s23_2, dtype=DTYPE))
c12 = np.cos(t12)
s12 = np.sin(t12)
@@ -98,9 +98,9 @@ def angles_to_u(bsm_angles):
c23 = np.cos(t23)
s23 = np.sin(t23)
- p1 = np.array([[1 , 0 , 0] , [0 , c23 , s23] , [0 , -s23 , c23]] , dtype=np.complex128)
- p2 = np.array([[c13 , 0 , s13*np.exp(-1j*dcp)] , [0 , 1 , 0] , [-s13*np.exp(1j*dcp) , 0 , c13]] , dtype=np.complex128)
- p3 = np.array([[c12 , s12 , 0] , [-s12 , c12 , 0] , [0 , 0 , 1]] , dtype=np.complex128)
+ p1 = np.array([[1 , 0 , 0] , [0 , c23 , s23] , [0 , -s23 , c23]] , dtype=CDTYPE)
+ p2 = np.array([[c13 , 0 , s13*np.exp(-1j*dcp)] , [0 , 1 , 0] , [-s13*np.exp(1j*dcp) , 0 , c13]] , dtype=CDTYPE)
+ p3 = np.array([[c12 , s12 , 0] , [-s12 , c12 , 0] , [0 , 0 , 1]] , dtype=CDTYPE)
u = np.dot(np.dot(p1, p2), p3)
return u
@@ -142,15 +142,15 @@ def cardano_eqn(ham):
a = -np.trace(ham)
b = (0.5) * ((np.trace(ham))**2 - np.trace(np.dot(ham, ham)))
- c = -linalg.det(ham)
+ c = -np.linalg.det(ham)
Q = (1/9.) * (a**2 - 3*b)
R = (1/54.) * (2*a**3 - 9*a*b + 27*c)
theta = np.arccos(R / np.sqrt(Q**3))
E1 = -2 * np.sqrt(Q) * np.cos(theta/3.) - (1/3.)*a
- E2 = -2 * np.sqrt(Q) * np.cos((theta - 2.*np.pi)/3.) - (1/3.)*a
- E3 = -2 * np.sqrt(Q) * np.cos((theta + 2.*np.pi)/3.) - (1/3.)*a
+ E2 = -2 * np.sqrt(Q) * np.cos((theta - 2.*PI)/3.) - (1/3.)*a
+ E3 = -2 * np.sqrt(Q) * np.cos((theta + 2.*PI)/3.) - (1/3.)*a
A1 = ham[1][2] * (ham[0][0] - E1) - ham[1][0]*ham[0][2]
A2 = ham[1][2] * (ham[0][0] - E2) - ham[1][0]*ham[0][2]
@@ -164,9 +164,9 @@ def cardano_eqn(ham):
C2 = ham[1][0] * (ham[2][2] - E2) - ham[1][2]*ham[2][0]
C3 = ham[1][0] * (ham[2][2] - E3) - ham[1][2]*ham[2][0]
- N1 = np.sqrt(abs(A1*B1)**2 + abs(A1*C1)**2 + abs(B1*C1)**2)
- N2 = np.sqrt(abs(A2*B2)**2 + abs(A2*C2)**2 + abs(B2*C2)**2)
- N3 = np.sqrt(abs(A3*B3)**2 + abs(A3*C3)**2 + abs(B3*C3)**2)
+ N1 = np.sqrt(np.abs(A1*B1, dtype=DTYPE)**2 + np.abs(A1*C1, dtype=DTYPE)**2 + np.abs(B1*C1, dtype=DTYPE)**2)
+ N2 = np.sqrt(np.abs(A2*B2, dtype=DTYPE)**2 + np.abs(A2*C2, dtype=DTYPE)**2 + np.abs(B2*C2, dtype=DTYPE)**2)
+ N3 = np.sqrt(np.abs(A3*B3, dtype=DTYPE)**2 + np.abs(A3*C3, dtype=DTYPE)**2 + np.abs(B3*C3, dtype=DTYPE)**2)
mm = np.array([
[np.conjugate(B1)*C1 / N1, np.conjugate(B2)*C2 / N2, np.conjugate(B3)*C3 / N3],
@@ -195,7 +195,7 @@ def normalise_fr(fr):
array([ 0.16666667, 0.33333333, 0.5 ])
"""
- return np.array(fr) / float(np.sum(fr))
+ return np.array(fr) / np.sum(fr)
def estimate_scale(args):
@@ -300,7 +300,7 @@ NUFIT_U = angles_to_u((0.307, (1-0.02195)**2, 0.565, 3.97935))
def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
sm_u=NUFIT_U, no_bsm=False, fix_mixing=False,
fix_mixing_almost=False, fix_scale=False, scale=None,
- check_uni=True, epsilon=1e-9):
+ check_uni=True, epsilon=1e-7):
"""Construct the BSM mixing matrix from the BSM parameters.
Parameters
@@ -393,8 +393,8 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
if check_uni:
tu = test_unitarity(eg_vector)
- if not abs(np.trace(tu) - 3.) < epsilon or \
- not abs(np.sum(tu) - 3.) < epsilon:
+ if not np.abs(np.trace(tu) - 3.) < epsilon or \
+ not np.abs(np.sum(tu) - 3.) < epsilon:
raise AssertionError(
'Matrix is not unitary!\neg_vector\n{0}\ntest '
'u\n{1}'.format(eg_vector, tu)
@@ -427,7 +427,7 @@ def test_unitarity(x, prnt=False):
[ 0., 0., 1.]])
"""
- f = abs(np.dot(x, x.conj().T))
+ f = np.abs(np.dot(x, x.conj().T), dtype=DTYPE)
if prnt:
print 'Unitarity test:\n{0}'.format(f)
return f
@@ -456,7 +456,8 @@ def u_to_fr(source_fr, matrix):
"""
composition = np.einsum(
- 'ai, bi, a -> b', abs(matrix)**2, abs(matrix)**2, source_fr
+ 'ai, bi, a -> b',
+ np.abs(matrix)**2, np.abs(matrix)**2, source_fr,
)
ratio = composition / np.sum(source_fr)
return ratio
diff --git a/utils/gf.py b/utils/gf.py
index 59d1033..17ac029 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -16,6 +16,7 @@ import GolemFitPy as gf
from utils.enums import DataType, SteeringCateg
from utils.misc import enum_parse, thread_factors
+from utils.param import ParamSet
def fit_flags(llh_paramset):
@@ -23,9 +24,6 @@ def fit_flags(llh_paramset):
# False means it's not fixed in minimization
'astroFlavorAngle1' : True,
'astroFlavorAngle2' : True,
- 'astroENorm' : True,
- 'astroMuNorm' : True,
- 'astroTauNorm' : True,
'convNorm' : True,
'promptNorm' : True,
'muonNorm' : True,
@@ -44,9 +42,14 @@ def fit_flags(llh_paramset):
'astroDeltaGammaSec' : True
}
flags = gf.FitParametersFlag()
+ gf_nuisance = []
for param in llh_paramset:
- flags.__setattr__(param.name, False)
- return flags
+ if param.name in default_flags:
+ print 'Setting param {0:<15} to float in the ' \
+ 'minimisation'.format(param.name)
+ flags.__setattr__(param.name, False)
+ gf_nuisance.append(param)
+ return flags, ParamSet(gf_nuisance)
def steering_params(args):
@@ -68,7 +71,7 @@ def set_up_as(fitter, params):
print 'Injecting the model', params
asimov_params = gf.FitParameters(gf.sampleTag.HESE)
for parm in params:
- asimov_params.__setattr__(parm.name, parm.value)
+ asimov_params.__setattr__(parm.name, float(parm.value))
fitter.SetupAsimov(asimov_params)
@@ -84,7 +87,7 @@ def setup_fitter(args, asimov_paramset):
def get_llh(fitter, params):
fitparams = gf.FitParameters(gf.sampleTag.HESE)
for parm in params:
- fitparams.__setattr__(parm.name, parm.value)
+ fitparams.__setattr__(parm.name, float(parm.value))
llh = -fitter.EvalLLH(fitparams)
return llh
@@ -93,9 +96,10 @@ def get_llh_freq(fitter, params):
print 'setting to {0}'.format(params)
fitparams = gf.FitParameters(gf.sampleTag.HESE)
for parm in params:
- fitparams.__setattr__(parm.name, parm.value)
+ fitparams.__setattr__(parm.name, float(parm.value))
fitter.SetFitParametersSeed(fitparams)
- return -fitter.MinLLH().likelihood
+ llh = -fitter.MinLLH().likelihood
+ return llh
def data_distributions(fitter):
diff --git a/utils/likelihood.py b/utils/likelihood.py
index 04828e8..6387a1e 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -12,7 +12,7 @@ from __future__ import absolute_import, division
from functools import partial
import numpy as np
-from scipy.stats import multivariate_normal
+from scipy.stats import multivariate_normal, rv_continuous
import GolemFitPy as gf
@@ -22,15 +22,16 @@ from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg
from utils.misc import enum_parse
-def multi_gaussian(fr, fr_bf, sigma):
- """Multivariate gaussian likelihood."""
- cov_fr = np.identity(3) * sigma
- return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr))
+class Gaussian(rv_continuous):
+ """Gaussian for one dimension."""
+ def _pdf(self, x, mu, sig):
+ return (1./np.sqrt(2*np.pi*sig**2))*np.exp(-((x-mu)**2)/(2*sig**2))
-def log_gaussian(x, mu, sig):
- """Log gaussian for one dimension."""
- return np.log((1./np.sqrt(2*np.pi*sig**2))*np.exp(-((x-mu)**2)/(2*sig**2)))
+def multi_gaussian(fr, fr_bf, sigma, offset=-320):
+ """Multivariate gaussian likelihood."""
+ cov_fr = np.identity(3) * sigma
+ return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr)) + offset
def likelihood_argparse(parser):
@@ -46,6 +47,11 @@ def likelihood_argparse(parser):
def lnprior(theta, paramset):
"""Priors on theta."""
+ if len(theta) != len(paramset):
+ raise AssertionError(
+ 'Length of MCMC scan is not the same as the input '
+ 'params\ntheta={0}\nparamset={1}'.format(theta, paramset)
+ )
for idx, param in enumerate(paramset):
param.value = theta[idx]
ranges = paramset.ranges
@@ -57,9 +63,13 @@ def lnprior(theta, paramset):
prior = 0
for param in paramset:
if param.prior is PriorsCateg.GAUSSIAN:
- prior += log_gaussian(
+ prior += Gaussian().logpdf(
param.nominal_value, param.value, param.std
)
+ elif param.prior is PriorsCateg.HALFGAUSS:
+ prior += Gaussian().logpdf(
+ param.nominal_value, param.value, param.std
+ ) + Gaussian().logcdf(1, param.value, param.std)
return prior
@@ -68,7 +78,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
if len(theta) != len(llh_paramset):
raise AssertionError(
'Length of MCMC scan is not the same as the input '
- 'params\ntheta={0}\nmcmc_paramset]{1}'.format(theta, llh_paramset)
+ 'params\ntheta={0}\nparamset]{1}'.format(theta, llh_paramset)
)
for idx, param in enumerate(llh_paramset):
param.value = theta[idx]
diff --git a/utils/misc.py b/utils/misc.py
index 59c5edb..970c693 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -22,6 +22,11 @@ import numpy as np
from utils.enums import Likelihood
+DTYPE = np.float128
+CDTYPE = np.complex128
+PI = np.arccos(DTYPE(-1))
+
+
class SortingHelpFormatter(argparse.HelpFormatter):
"""Sort argparse help options alphabetically."""
def add_arguments(self, actions):
diff --git a/utils/multinest.py b/utils/multinest.py
index 426a951..005a43a 100644
--- a/utils/multinest.py
+++ b/utils/multinest.py
@@ -16,26 +16,11 @@ import numpy as np
from pymultinest import analyse, run
from utils import likelihood
-from utils.misc import make_dir, parse_bool
+from utils.misc import gen_outfile_name, make_dir, parse_bool
-def CubePrior(cube, ndim, n_params, mn_paramset):
- if ndim != len(mn_paramset):
- raise AssertionError(
- 'Length of MultiNest scan paramset is not the same as the input '
- 'params\ncube={0}\nmn_paramset]{1}'.format(cube, mn_paramset)
- )
- pranges = mn_paramset.seeds
- for i in xrange(ndim):
- mn_paramset[i].value = (pranges[i][1]-pranges[i][0])*cube[i] + pranges[i][0]
- prior = 0
- for parm in mn_paramset:
- if parm.prior is PriorsCateg.GAUSSIAN:
- prior_penatly += multivariate_normal(
- parm.nominal_value, mean=parm.value, cov=parm.std
- )
- print 'prior', prior
- return prior
+def CubePrior(cube, ndim, n_params):
+ pass
def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
@@ -63,18 +48,6 @@ def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
def mn_argparse(parser):
parser.add_argument(
- '--mn-run', type=parse_bool, default='True',
- help='Run MultiNest'
- )
- parser.add_argument(
- '--mn-bins', type=int, default=10,
- help='Number of bins for the Bayes factor plot'
- )
- parser.add_argument(
- '--mn-eval-bin', type=str, default='all',
- help='Which bin to evalaute for Bayes factor plot'
- )
- parser.add_argument(
'--mn-live-points', type=int, default=3000,
help='Number of live points for MultiNest evaluations'
)
@@ -104,8 +77,8 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
fitter = fitter
)
- prefix = './mnrun/DIM{0}/{1:>019}_'.format(
- args.dimension, np.random.randint(0, 2**63)
+ prefix = './mnrun/DIM{0}/{1}_{2:>019}_'.format(
+ args.dimension, gen_outfile_name(args), np.random.randint(0, 2**63)
)
make_dir(prefix)
print 'Running evidence calculation for {0}'.format(prefix)
diff --git a/utils/param.py b/utils/param.py
index 13f1a63..4d8106e 100644
--- a/utils/param.py
+++ b/utils/param.py
@@ -205,6 +205,13 @@ class ParamSet(Sequence):
else:
return ParamSet([io[1] for io in ps])
+ def remove_params(self, params):
+ rm_paramset = []
+ for parm in self.params:
+ if parm.name not in params.names:
+ rm_paramset.append(parm)
+ return ParamSet(rm_paramset)
+
def get_paramsets(args, nuisance_paramset):
"""Make the paramsets for generating the Asmimov MC sample and also running
@@ -216,7 +223,7 @@ def get_paramsets(args, nuisance_paramset):
llh_paramset.extend(
[x for x in nuisance_paramset.from_tag(ParamTag.SM_ANGLES)]
)
- if args.likelihood is Likelihood.GOLEMFIT:
+ if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
gf_nuisance = [x for x in nuisance_paramset.from_tag(ParamTag.NUISANCE)]
asimov_paramset.extend(gf_nuisance)
llh_paramset.extend(gf_nuisance)
diff --git a/utils/plot.py b/utils/plot.py
index 6ae8dc2..0c431cf 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -233,10 +233,11 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset):
def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
"""Make MultiNest factor or LLH value plot."""
- print "Making MultiNest Factor 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
if args.stat_method is StatCateg.BAYESIAN:
@@ -244,9 +245,9 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
null = statistic[min_idx]
reduced_ev = -(statistic - null)
elif args.stat_method is StatCateg.FREQUENTIST:
- min_idx = np.argmin(statistic)
+ min_idx = np.argmin(scales)
null = statistic[min_idx]
- reduced_ev = 2*(statistic - null)
+ reduced_ev = -2*(statistic - null)
fig = plt.figure(figsize=(7, 5))
ax = fig.add_subplot(111)
@@ -256,7 +257,7 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
if args.stat_method is StatCateg.BAYESIAN:
ax.set_ylabel(r'Bayes Factor')
elif args.stat_method is StatCateg.FREQUENTIST:
- ax.set_ylabel(r'$\Delta {\rm LLH}$')
+ ax.set_ylabel(r'$-2\Delta {\rm LLH}$')
ax.plot(scales, reduced_ev)
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-17 09:50:54 +0000