Tue 20 Nov 23:28:31 CST 2018

8 files changed, 62 insertions, 41 deletions

diff --git a/plot_llh/calc_fr_MCMC.py b/plot_llh/calc_fr_MCMC.py
index ad75043..1bf016a 100644
--- a/plot_llh/calc_fr_MCMC.py
+++ b/plot_llh/calc_fr_MCMC.py
@@ -13,8 +13,8 @@ from utils.enums import MixingScenario
 
 binning = np.logspace(np.log10(6e4), np.log10(1e7), 21)
 dimension = 6
-source = [1, 0, 0]
-scenario = MixingScenario.T23
+source = [0, 1, 0]
+scenario = MixingScenario.T13
 
 def get_fr(theta, source, binning, dimension, scenario):
     sm_mixings = theta[:6]
@@ -56,7 +56,7 @@ def get_fr(theta, source, binning, dimension, scenario):
         intergrated_measured_flux
     fr = averaged_measured_flux / np.sum(averaged_measured_flux)
 
-    return fr
+    return map(float, fr)
 
 if len(sys.argv)< 2:
     print sys.argv
diff --git a/sens.py b/sens.py
index a97370b..4775c21 100755
--- a/sens.py
+++ b/sens.py
@@ -38,7 +38,7 @@ def define_nuisance():
     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.950, 0.961],   ranges=[0., 1.],      std=0.00147, tex=r'c_{13}^4', 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(
@@ -108,6 +108,8 @@ def process_args(args):
        args.likelihood is Likelihood.GOLEMFIT:
         args.likelihood = Likelihood.GF_FREQ
 
+    args.burnin = False
+
 
 def parse_args(args=None):
     """Parse command line arguments"""
@@ -273,13 +275,17 @@ def main():
                     print '|||| SCALE = {0:.0E}'.format(np.power(10, sc))
                     scale.value = sc
                     if args.stat_method is StatCateg.BAYESIAN:
+                        identifier = 'b{0}_{1}_sce{2}_sca{3}_an{4}'.format(
+                            args.sens_eval_bin, args.sens_bins, idx_scen, sc, idx_an
+                        )
                         try:
                             stat = mn_utils.mn_evidence(
                                 mn_paramset     = sens_paramset,
                                 llh_paramset    = llh_paramset,
                                 asimov_paramset = asimov_paramset,
                                 args            = args,
-                                fitter          = fitter
+                                fitter          = fitter,
+                                identifier      = identifier
                             )
                         except:
                             print 'Failed run, continuing'
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index be99d8b..0a48baa 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -29,19 +29,20 @@ GLOBAL_PARAMS.update(dict(
     sens_run      = 'True',
     run_method    = 'fixed_angle', # full, fixed_angle, corr_angle
     stat_method   = 'bayesian',
-    sens_bins     = 20,
+    sens_bins     = 10,
     seed          = None
 ))
 
 # MultiNest
 GLOBAL_PARAMS.update(dict(
     mn_live_points = 1000,
-    mn_tolerance   = 0.1,
+    # mn_tolerance   = 0.1,
+    mn_tolerance   = 0.3,
     mn_output      = './mnrun'
 ))
 
 # FR
-# dimension         = [3]
+# dimension         = [6]
 # dimension         = [3, 6]
 dimension         = [3, 4, 5, 6, 7, 8]
 GLOBAL_PARAMS.update(dict(
@@ -78,6 +79,7 @@ outfile = 'dagman_FR_SENS_{0}_{1}_{2}_{3}'.format(
     GLOBAL_PARAMS['likelihood'], GLOBAL_PARAMS['data']
 )
 # outfile += '_seed2'
+# outfile += '_tol03'
 outfile += '.submit'
 golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
 condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
@@ -102,6 +104,8 @@ with open(outfile, 'w') as f:
             if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
                 output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
             # output += 'seed2/'
+            # output += 'mn_noverlap/'
+            # output += 'tol_03/'
             for r in xrange(sens_runs):
                 print 'run', r
                 f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
@@ -121,6 +125,7 @@ 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
+                # break
 
         # for frs in full_scan_mfr:
         #     print 'frs', frs
diff --git a/submitter/sens_submit.sub b/submitter/sens_submit.sub
index 5b34dfd..6c92337 100644
--- a/submitter/sens_submit.sub
+++ b/submitter/sens_submit.sub
@@ -1,5 +1,5 @@
 Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/sens.py
-Arguments = "--ast $(ast) --data $(data) --dimension $(dimension) --fix-mixing $(fix_mixing) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --no-bsm $(no_bsm) --outfile $(outfile) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning) --fix-mixing-almost $(fix_mixing_almost) --sens-run $(sens_run) --run-method $(run_method) --stat-method $(stat_method) --sens-bins $(sens_bins) --sens-eval-bin $(sens_eval_bin) --mn-live-points $(mn_live_points) --mn-tolerance $(mn_tolerance) --mn-output $(mn_output) --plot-statistic $(plot_statistic) --fold-index $(fold_index)"
+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) --fold-index $(fold_index) --save-measured-fr $(save_measured_fr) --output-measured-fr=$(output_measured_fr)"
 
 # All logs will go to a single file 
 log    = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
diff --git a/utils/gf.py b/utils/gf.py
index 6ce7863..92e64f3 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -19,7 +19,7 @@ except:
     print 'Running without GolemFit'
     pass
 
-from utils.enums import DataType, SteeringCateg
+from utils.enums import DataType, Likelihood, SteeringCateg
 from utils.misc import enum_parse, thread_factors
 from utils.param import ParamSet
 
@@ -62,10 +62,16 @@ def steering_params(args):
     params = gf.SteeringParams(gf.sampleTag.MagicTau)
     params.quiet = False
     params.fastmode = True
+    # params.fastmode = False
     params.simToLoad= steering_categ.name.lower()
     params.evalThreads = args.threads
     # params.evalThreads = thread_factors(args.threads)[1]
 
+    if args.likelihood is Likelihood.GOLEMFIT:
+        params.frequentist = false;
+    elif args.likelihood is Likelihood.GF_FREQ:
+        params.frequentist = true;
+
     # For Tianlu
     # params.years = [999]
 
diff --git a/utils/likelihood.py b/utils/likelihood.py
index 7079d52..f4404c4 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -72,8 +72,6 @@ def lnprior(theta, paramset):
             prior += Gaussian().logpdf(
                 param.nominal_value, param.value, param.std
             )
-            print 'prioring param', param.name, '=', param.value
-            print 'prior', prior
         elif param.prior is PriorsCateg.HALFGAUSS:
             prior += Gaussian().logpdf(
                 param.nominal_value, param.value, param.std
@@ -175,15 +173,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
             intergrated_measured_flux
         fr = averaged_measured_flux / np.sum(averaged_measured_flux)
 
-    if args.save_measured_fr and args.burnin is False:
-        n = gen_identifier(args) + '.txt'
-        with open(args.output_measured_fr + n, 'a') as f:
-            f.write(r'{0:.3f} {1:.3f} {2:.3f} {3:.1f}'.format(
-                float(fr[0]), float(fr[1]), float(fr[2]),
-                llh_paramset['logLam'].value
-            ))
-            f.write('\n')
-
     flavour_angles = fr_utils.fr_to_angles(fr)
     # print 'flavour_angles', map(float, flavour_angles)
     for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)):
@@ -198,6 +187,18 @@ 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:
         llh = gf_utils.get_llh_freq(fitter, hypo_paramset)
+
+    if args.save_measured_fr and args.burnin is False:
+        n = gen_identifier(args) + '.txt'
+        with open(args.output_measured_fr + n, 'a') as f:
+            f.write(r'{0:.3f} {1:.3f} {2:.3f}'.format(
+                float(fr[0]), float(fr[1]), float(fr[2])
+            ))
+            for p in llh_paramset:
+                f.write(r' {0:.3f}'.format(p.value))
+            f.write('     LLH = {0:.3f}'.format(llh))
+            f.write('\n')
+
     return llh
 
 
@@ -205,10 +206,6 @@ def ln_prob(theta, args, asimov_paramset, llh_paramset, fitter):
     lp = lnprior(theta, paramset=llh_paramset)
     if not np.isfinite(lp):
         return -np.inf
-    llh = triangle_llh(
-        theta, args=args, asimov_paramset=asimov_paramset,
-        llh_paramset=llh_paramset, fitter=fitter
-    )
     return lp + triangle_llh(
         theta, args=args, asimov_paramset=asimov_paramset,
         llh_paramset=llh_paramset, fitter=fitter
diff --git a/utils/multinest.py b/utils/multinest.py
index e9eece9..fdd87cd 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_identifier, make_dir
+from utils.misc import gen_identifier, make_dir, solve_ratio
 
 
 def CubePrior(cube, ndim, n_params):
@@ -63,7 +63,8 @@ def mn_argparse(parser):
     )
 
 
-def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
+def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter,
+                identifier='mn'):
     """Run the MultiNest algorithm to calculate the evidence."""
     n_params = len(mn_paramset)
 
@@ -79,8 +80,14 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
         fitter          = fitter
     )
 
-    prefix = './mnrun/DIM{0}/{1}_{2}_{3:>010}_'.format(
-        args.dimension, args.likelihood, gen_identifier(args), np.random.randint(0, 2**30)
+    # prefix = './mnrun/DIM{0}/{1}_{2}_{3:>010}_'.format(
+    #     args.dimension, args.likelihood, gen_identifier(args), np.random.randint(0, 2**30)
+    # )
+    llh = '{0}'.format(args.likelihood).split('.')[1]
+    data = '{0}'.format(args.data).split('.')[1]
+    sr1, sr2, sr3 = solve_ratio(args.source_ratio)
+    prefix = './mnrun/DIM{0}/{1}/{2}/s{3}{4}{5}/{6}'.format(
+        args.dimension, data, llh, sr1, sr2, sr3, identifier
     )
     make_dir(prefix)
     print 'Running evidence calculation for {0}'.format(prefix)
diff --git a/utils/plot.py b/utils/plot.py
index 53a4a84..6cbc4f5 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -627,18 +627,18 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
     plt.setp(legend.get_title(), fontsize='11')
     legend.get_frame().set_linestyle('-')
 
-    # if show_data: ybound = 0.595
-    # else: ybound = 0.73
-    # if args.data is DataType.REAL and show_data:
-    #     # fig.text(0.295, 0.684, 'IceCube Preliminary', color='red', fontsize=13,
-    #     fig.text(0.278, ybound, r'\bf IceCube Preliminary', color='red', fontsize=13,
-    #              ha='center', va='center', zorder=11)
-    # elif args.data is DataType.REALISATION:
-    #     fig.text(0.278, ybound-0.05, r'\bf IceCube Simulation', color='red', fontsize=13,
-    #              ha='center', va='center', zorder=11)
-    # else:
-    #     fig.text(0.278, ybound, r'\bf IceCube Simulation', color='red', fontsize=13,
-    #              ha='center', va='center', zorder=11)
+    if show_data: ybound = 0.595
+    else: ybound = 0.73
+    if args.data is DataType.REAL and show_data:
+        # fig.text(0.295, 0.684, 'IceCube Preliminary', color='red', fontsize=13,
+        fig.text(0.278, ybound, r'\bf IceCube Preliminary', color='red', fontsize=13,
+                 ha='center', va='center', zorder=11)
+    elif args.data is DataType.REALISATION:
+        fig.text(0.278, ybound-0.05, r'\bf IceCube Simulation', color='red', fontsize=13,
+                 ha='center', va='center', zorder=11)
+    else:
+        fig.text(0.278, ybound, r'\bf IceCube Simulation', color='red', fontsize=13,
+                 ha='center', va='center', zorder=11)
 
     misc_utils.make_dir(outfile)
     for of in outformat: