Mon Apr 16 16:57:08 CDT 2018

8 files changed, 417 insertions, 62 deletions

diff --git a/.gitignore b/.gitignore
index 8728cfe..7dd464f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,3 +8,4 @@ dagman_FR*
 submitter/logs/
 mnrun_*
 *.png
+mnrun/
diff --git a/fr.py b/fr.py
index 90c2e94..906c0bf 100755
--- a/fr.py
+++ b/fr.py
@@ -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')