Sun Apr 22 23:18:44 CDT 2018

9 files changed, 541 insertions, 309 deletions

diff --git a/utils/enums.py b/utils/enums.py
index 90ca17d..359e104 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -47,6 +47,17 @@ class MCMCSeedType(Enum):
     GAUSSIAN = 2
 
 
+class SensitivityCateg(Enum):
+    FULL        = 1
+    FIXED_ANGLE = 2
+    CORR_ANGLE  = 3
+
+
+class StatCateg(Enum):
+    BAYESIAN    = 1
+    FREQUENTIST = 2
+
+
 class SteeringCateg(Enum):
     P2_0 = 1
     P2_1 = 2
diff --git a/utils/fr.py b/utils/fr.py
index 8f71f78..6a405a7 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -9,11 +9,15 @@ Useful functions for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import sys
+import argparse
+from functools import partial
 
 import numpy as np
 from scipy import linalg
 
+from utils.enums import EnergyDependance
+from utils.misc import enum_parse, parse_bool
+
 
 MASS_EIGENVALUES = [7.40E-23, 2.515E-21]
 """SM mass eigenvalues"""
@@ -194,6 +198,83 @@ def normalise_fr(fr):
     return np.array(fr) / float(np.sum(fr))
 
 
+def estimate_scale(args, mass_eigenvalues=MASS_EIGENVALUES):
+    """Estimate the scale at which new physics will enter."""
+    if args.energy_dependance is EnergyDependance.MONO:
+        scale = np.power(
+            10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \
+            np.log10(args.energy**(args.dimension-3))
+        )
+    elif args.energy_dependance is EnergyDependance.SPECTRAL:
+        scale = np.power(
+            10, np.round(
+                np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \
+                - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3))
+            )
+        )
+    return scale
+
+
+def fr_argparse(parser):
+    parser.add_argument(
+        '--measured-ratio', type=int, nargs=3, default=[1, 1, 1],
+        help='Set the central value for the measured flavour ratio at IceCube'
+    )
+    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'
+    )
+    parser.add_argument(
+        '--no-bsm', type=parse_bool, default='False',
+        help='Turn off BSM terms'
+    )
+    parser.add_argument(
+        '--dimension', type=int, default=3,
+        help='Set the new physics dimension to consider'
+    )
+    parser.add_argument(
+        '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance),
+        choices=EnergyDependance,
+        help='Type of energy dependance to use in the BSM fit'
+    )
+    parser.add_argument(
+        '--spectral-index', default=-2, type=int,
+        help='Spectral index for spectral energy dependance'
+    )
+    parser.add_argument(
+        '--binning', default=[1e4, 1e7, 5], type=float, nargs=3,
+        help='Binning for spectral energy dependance'
+    )
+    parser.add_argument(
+        '--fix-source-ratio', type=parse_bool, default='False',
+        help='Fix the source flavour ratio'
+    )
+    parser.add_argument(
+        '--fix-mixing', type=parse_bool, default='False',
+        help='Fix all mixing parameters to bi-maximal mixing'
+    )
+    parser.add_argument(
+        '--fix-mixing-almost', type=parse_bool, default='False',
+        help='Fix all mixing parameters except s23'
+    )
+    parser.add_argument(
+        '--fix-scale', type=parse_bool, default='False',
+        help='Fix the new physics scale'
+    )
+    parser.add_argument(
+        '--scale', type=float, default=1e-30,
+        help='Set the new physics scale'
+    )
+    parser.add_argument(
+        '--scale-region', type=float, default=1e10,
+        help='Set the size of the box to scan for the scale'
+    )
+    parser.add_argument(
+        '--energy', type=float, default=1000,
+        help='Set the energy scale'
+    )
+
+
 def fr_to_angles(ratios):
     """Convert from flavour ratio into the angular projection of the flavour
     ratios.
@@ -218,7 +299,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,
                    nufit_u=NUFIT_U, no_bsm=False, fix_mixing=False,
                    fix_mixing_almost=False, fix_scale=False, scale=None,
-                   check_uni=True):
+                   check_uni=True, epsilon=1e-9):
     """Construct the BSM mixing matrix from the BSM parameters.
 
     Parameters
@@ -311,8 +392,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.) < 1e-5 or \
-           not abs(np.sum(tu) - 3.) < 1e-5:
+        if not abs(np.trace(tu) - 3.) < epsilon or \
+           not abs(np.sum(tu) - 3.) < epsilon:
             raise AssertionError(
                 'Matrix is not unitary!\neg_vector\n{0}\ntest '
                 'u\n{1}'.format(eg_vector, tu)
@@ -378,4 +459,3 @@ def u_to_fr(source_fr, matrix):
     )
     ratio = composition / np.sum(source_fr)
     return ratio
-
diff --git a/utils/gf.py b/utils/gf.py
index 3f770eb..20afc75 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -9,7 +9,6 @@ Useful GolemFit wrappers for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import argparse
 import socket
 from functools import partial
 
diff --git a/utils/likelihood.py b/utils/likelihood.py
index c1208ab..1981c72 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -9,7 +9,6 @@ Likelihood functions for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import argparse
 from functools import partial
 
 import numpy as np
@@ -34,6 +33,10 @@ def likelihood_argparse(parser):
         '--likelihood', default='gaussian', type=partial(enum_parse, c=Likelihood),
         choices=Likelihood, help='likelihood contour'
     )
+    parser.add_argument(
+        '--sigma-ratio', type=float, default=0.01,
+        help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH'
+    )
 
 
 def lnprior(theta, paramset):
@@ -46,17 +49,17 @@ def lnprior(theta, paramset):
     return 0.
 
 
-def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
+def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
     """-Log likelihood function for a given theta."""
-    if len(theta) != len(mcmc_paramset):
+    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, mcmc_paramset)
+            'params\ntheta={0}\nmcmc_paramset]{1}'.format(theta, llh_paramset)
         )
-    for idx, param in enumerate(mcmc_paramset):
+    for idx, param in enumerate(llh_paramset):
         param.value = theta[idx]
     hypo_paramset = asimov_paramset
-    for param in mcmc_paramset.from_tag(ParamTag.NUISANCE):
+    for param in llh_paramset.from_tag(ParamTag.NUISANCE):
         hypo_paramset[param.name].value = param.value
 
     if args.energy_dependance is EnergyDependance.SPECTRAL:
@@ -67,14 +70,14 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
         if args.energy_dependance is EnergyDependance.MONO:
             source_flux = args.source_ratio
         elif args.energy_dependance is EnergyDependance.SPECTRAL:
-	    source_flux = np.array(
-		[fr * np.power(bin_centers, args.spectral_index)
-		 for fr in args.source_ratio]
-	    ).T
+            source_flux = np.array(
+                [fr * np.power(bin_centers, args.spectral_index)
+                 for fr in args.source_ratio]
+            ).T
     else:
         if args.energy_dependance is EnergyDependance.MONO:
             source_flux = fr_utils.angles_to_fr(
-                mcmc_paramset.from_tag(ParamTag.SRCANGLES, values=True)
+                llh_paramset.from_tag(ParamTag.SRCANGLES, values=True)
             )
         elif args.energy_dependance is EnergyDependance.SPECTRAL:
             source_flux = np.array(
@@ -82,7 +85,7 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
                  for fr in fr_utils.angles_to_fr(theta[-2:])]
             ).T
 
-    bsm_angles = mcmc_paramset.from_tag(
+    bsm_angles = llh_paramset.from_tag(
         [ParamTag.SCALE, ParamTag.MMANGLES], values=True
     )
 
@@ -134,11 +137,11 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
         return gf_utils.get_llh_freq(fitter, hypo_paramset)
 
 
-def ln_prob(theta, args, fitter, asimov_paramset, mcmc_paramset):
-    lp = lnprior(theta, paramset=mcmc_paramset)
+def ln_prob(theta, args, fitter, asimov_paramset, llh_paramset):
+    lp = lnprior(theta, paramset=llh_paramset)
     if not np.isfinite(lp):
         return -np.inf
     return lp + triangle_llh(
         theta, args=args, asimov_paramset=asimov_paramset,
-        mcmc_paramset=mcmc_paramset, fitter=fitter
+        llh_paramset=llh_paramset, fitter=fitter
     )
diff --git a/utils/mcmc.py b/utils/mcmc.py
index c712c3a..97b77dd 100644
--- a/utils/mcmc.py
+++ b/utils/mcmc.py
@@ -9,7 +9,6 @@ Useful functions to use an MCMC for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import argparse
 from functools import partial
 
 import emcee
@@ -72,6 +71,14 @@ def mcmc_argparse(parser):
         type=partial(enum_parse, c=MCMCSeedType), choices=MCMCSeedType,
         help='Type of distrbution to make the initial MCMC seed'
     )
+    parser.add_argument(
+        '--plot-angles', type=misc_utils.parse_bool, default='True',
+        help='Plot MCMC triangle in the angles space'
+    )
+    parser.add_argument(
+        '--plot-elements', type=misc_utils.parse_bool, default='False',
+        help='Plot MCMC triangle in the mixing elements space'
+    )
 
 
 def flat_seed(paramset, nwalkers):
diff --git a/utils/misc.py b/utils/misc.py
index f0c1ad4..59c5edb 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -9,7 +9,7 @@ Misc functions for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import os, sys
+import os
 import errno
 import multiprocessing
 
@@ -19,174 +19,7 @@ from operator import attrgetter
 
 import numpy as np
 
-from utils.enums import Likelihood, ParamTag
-
-
-class Param(object):
-    """Parameter class to store parameters.
-    """
-    def __init__(self, name, value, ranges, seed=None, std=None, tex=None, tag=None):
-        self._seed = None
-        self._ranges = None
-        self._tex = None
-        self._tag = None
-
-        self.name = name
-        self.value = value
-        self.seed = seed
-        self.ranges = ranges
-        self.std = std
-        self.tex = tex
-        self.tag = tag
-
-    @property
-    def ranges(self):
-        return tuple(self._ranges)
-
-    @ranges.setter
-    def ranges(self, values):
-        self._ranges = [val for val in values]
-
-    @property
-    def seed(self):
-        if self._seed is None: return self.ranges
-        return tuple(self._seed)
-
-    @seed.setter
-    def seed(self, values):
-        if values is None: return
-        self._seed = [val for val in values]
-
-    @property
-    def tex(self):
-        return r'{0}'.format(self._tex)
-
-    @tex.setter
-    def tex(self, t):
-        self._tex = t if t is not None else r'{\rm %s}' % self.name
-
-    @property
-    def tag(self):
-        return self._tag
-
-    @tag.setter
-    def tag(self, t):
-        if t is None: self._tag = ParamTag.NONE
-        else:
-            assert t in ParamTag
-            self._tag = t
-
-
-class ParamSet(Sequence):
-    """Container class for a set of parameters.
-    """
-    def __init__(self, *args):
-        param_sequence = []
-        for arg in args:
-            try:
-                param_sequence.extend(arg)
-            except TypeError:
-                param_sequence.append(arg)
-
-        if len(param_sequence) != 0:
-            # Disallow duplicated params
-            all_names = [p.name for p in param_sequence]
-            unique_names = set(all_names)
-            if len(unique_names) != len(all_names):
-                duplicates = set([x for x in all_names if all_names.count(x) > 1])
-                raise ValueError('Duplicate definitions found for param(s): ' +
-                                 ', '.join(str(e) for e in duplicates))
-
-        # Elements of list must be Param type
-        assert all([isinstance(x, Param) for x in param_sequence]), \
-                'All params must be of type "Param"'
-
-        self._params = param_sequence
-
-    def __len__(self):
-        return len(self._params)
-
-    def __getitem__(self, i):
-        if isinstance(i, int):
-            return self._params[i]
-        elif isinstance(i, basestring):
-            return self._by_name[i]
-
-    def __getattr__(self, attr):
-        try:
-            return super(ParamSet, self).__getattribute__(attr)
-        except AttributeError:
-            t, v, tb = sys.exc_info()
-            try:
-                return self[attr]
-            except KeyError:
-                raise t, v, tb
-
-    def __iter__(self):
-        return iter(self._params)
-
-    def __str__(self):
-        o = '\n'
-        for obj in self._params:
-            o += '== {0:<15} = {1:<15}, tag={2:<15}\n'.format(
-                obj.name, obj.value, obj.tag
-            )
-        return o
-
-    @property
-    def _by_name(self):
-        return {obj.name: obj for obj in self._params}
-
-    @property
-    def names(self):
-        return tuple([obj.name for obj in self._params])
-
-    @property
-    def labels(self):
-        return tuple([obj.tex for obj in self._params])
-
-    @property
-    def values(self):
-        return tuple([obj.value for obj in self._params])
-
-    @property
-    def seeds(self):
-        return tuple([obj.seed for obj in self._params])
-
-    @property
-    def ranges(self):
-        return tuple([obj.ranges for obj in self._params])
-
-    @property
-    def stds(self):
-        return tuple([obj.std for obj in self._params])
-
-    @property
-    def tags(self):
-        return tuple([obj.tag for obj in self._params])
-
-    @property
-    def params(self):
-        return self._params
-
-    def to_dict(self):
-        return {obj.name: obj.value for obj in self._params}
-
-    def from_tag(self, tag, values=False, index=False, invert=False):
-        if values and index: assert 0
-        tag = np.atleast_1d(tag)
-        if not invert:
-            ps = [(idx, obj) for idx, obj in enumerate(self._params)
-                  if obj.tag in tag]
-        else:
-            ps = [(idx, obj) for idx, obj in enumerate(self._params)
-                  if obj.tag not in tag]
-        if values:
-            return tuple([io[1].value for io in ps])
-        elif index:
-            return tuple([io[0] for io in ps])
-        else:
-            return ParamSet([io[1] for io in ps])
+from utils.enums import Likelihood
 
 
 class SortingHelpFormatter(argparse.HelpFormatter):
@@ -197,54 +30,32 @@ class SortingHelpFormatter(argparse.HelpFormatter):
 
 
 def gen_identifier(args):
-    mr = args.measured_ratio
-    si = args.sigma_ratio
+    f = '_DIM{0}'.format(args.dimension)
+    mr1, mr2, mr3 = args.measured_ratio
     if args.fix_source_ratio:
-        sr = args.source_ratio
+        sr1, sr2, sr3 = args.source_ratio
+        f += '_sfr_{0:03d}_{1:03d}_{2:03d}_mfr_{3:03d}_{4:03d}_{5:03d}'.format(
+            int(sr1*100), int(sr2*100), int(sr3*100),
+            int(mr1*100), int(mr2*100), int(mr3*100)
+        )
         if args.fix_mixing:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fix_mixing'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
-                int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension
-            )
+            f += '_fix_mixing'
         elif args.fix_mixing_almost:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fix_mixing_almost'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
-                int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension
-            )
+            f += '_fix_mixing_almost'
         elif args.fix_scale:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fixed_scale_{8}'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
-                int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension,
-                args.scale
-            )
-        else:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
-                int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension
-            )
+            f += '_fix_scale_{0}'.format(args.scale)
     else:
-        if args.fix_mixing:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fix_mixing'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
-                int(si*1000), args.dimension
-            )
-        elif args.fix_mixing_almost:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fix_mixing_almost'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
-                int(si*1000), args.dimension
-            )
-        elif args.fix_scale:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fixed_scale_{5}'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
-                int(si*1000), args.dimension, args.scale
-            )
-        else:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
-                int(si*1000), args.dimension
-            )
-    if args.likelihood is Likelihood.FLAT: out += '_flat'
-    return out
+        f += '_mfr_{3:03d}_{4:03d}_{5:03d}'.format(mr1, mr2, mr3)
+    if args.fix_mixing:
+        f += '_fix_mixing'
+    elif args.fix_mixing_almost:
+        f += '_fix_mixing_almost'
+    elif args.fix_scale:
+        f += '_fix_scale_{0}'.format(args.scale)
+    if args.likelihood is Likelihood.FLAT: f += '_flat'
+    elif args.likelihood is Likelihood.GAUSSIAN:
+        f += '_sigma_{0:03d}'.format(int(args.sigma_ratio*1000))
+    return f
 
 
 def gen_outfile_name(args):
@@ -314,6 +125,13 @@ def make_dir(outfile):
             raise
 
 
+def seed_parse(s):
+    if s.lower() == 'none':
+        return None
+    else:
+        return int(s)
+
+
 def thread_type(t):
     if t.lower() == 'max':
         return multiprocessing.cpu_count()
diff --git a/utils/multinest.py b/utils/multinest.py
new file mode 100644
index 0000000..3a7ab2d
--- /dev/null
+++ b/utils/multinest.py
@@ -0,0 +1,117 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : April 19, 2018
+
+"""
+Useful functions to use MultiNest for the BSM flavour ratio analysis
+"""
+
+from __future__ import absolute_import, division
+
+from functools import partial
+
+import numpy as np
+
+from pymultinest import analyse, run
+
+from utils.likelihood import triangle_llh
+from utils.misc import make_dir, parse_bool
+
+
+def CubePrior(cube, ndim, n_params):
+    # default are uniform priors
+    return ;
+
+
+def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
+           args, fitter):
+    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]
+    for pm in mn_paramset.names:
+        llh_paramset[pm].value = mn_paramset[pm].value
+    theta = llh_paramset.values
+    # print 'llh_paramset', llh_paramset
+    return triangle_llh(
+        theta=theta,
+        args=args,
+        asimov_paramset=asimov_paramset,
+        llh_paramset=llh_paramset,
+        fitter=fitter
+    )
+
+
+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'
+    )
+    parser.add_argument(
+        '--mn-tolerance', type=float, default=0.01,
+        help='Tolerance for MultiNest'
+    )
+    parser.add_argument(
+        '--mn-output', type=str, default='./mnrun/',
+        help='Folder to store MultiNest evaluations'
+    )
+    parser.add_argument(
+        '--plot-multinest', type=parse_bool, default='False',
+        help='Plot MultiNest evidence'
+    )
+
+
+def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
+    """Run the MultiNest algorithm to calculate the evidence."""
+    n_params = len(mn_paramset)
+
+    for n in mn_paramset.names:
+        llh_paramset[n].value = mn_paramset[n].value
+
+    lnProbEval = partial(
+        lnProb,
+        mn_paramset     = mn_paramset,
+        llh_paramset    = llh_paramset,
+        asimov_paramset = asimov_paramset,
+        args            = args,
+        fitter          = fitter
+    )
+
+    prefix = './mnrun/DIM{0}/{1:>019}_'.format(
+        args.dimension, np.random.randint(0, 2**63)
+    )
+    make_dir(prefix)
+    print 'Running evidence calculation for {0}'.format(prefix)
+    result = run(
+        LogLikelihood              = lnProbEval,
+        Prior                      = CubePrior,
+        n_dims                     = n_params,
+        n_live_points              = args.mn_live_points,
+        evidence_tolerance         = args.mn_tolerance,
+        outputfiles_basename       = prefix,
+        importance_nested_sampling = True,
+        resume                     = False,
+        verbose                    = True
+    )
+
+    analyser = analyse.Analyzer(
+        outputfiles_basename=prefix, n_params=n_params
+    )
+    return analyser.get_stats()['global evidence']
diff --git a/utils/param.py b/utils/param.py
new file mode 100644
index 0000000..bf230df
--- /dev/null
+++ b/utils/param.py
@@ -0,0 +1,235 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : April 19, 2018
+
+"""
+Param class and functions for the BSM flavour ratio analysis
+"""
+
+from __future__ import absolute_import, division
+
+import sys
+
+from collections import Sequence
+
+import numpy as np
+
+from utils.plot import get_units
+from utils.fr import fr_to_angles
+from utils.enums import Likelihood, ParamTag
+
+
+class Param(object):
+    """Parameter class to store parameters."""
+    def __init__(self, name, value, ranges, seed=None, std=None, tex=None, tag=None):
+        self._seed = None
+        self._ranges = None
+        self._tex = None
+        self._tag = None
+
+        self.name = name
+        self.value = value
+        self.seed = seed
+        self.ranges = ranges
+        self.std = std
+        self.tex = tex
+        self.tag = tag
+
+    @property
+    def ranges(self):
+        return tuple(self._ranges)
+
+    @ranges.setter
+    def ranges(self, values):
+        self._ranges = [val for val in values]
+
+    @property
+    def seed(self):
+        if self._seed is None: return self.ranges
+        return tuple(self._seed)
+
+    @seed.setter
+    def seed(self, values):
+        if values is None: return
+        self._seed = [val for val in values]
+
+    @property
+    def tex(self):
+        return r'{0}'.format(self._tex)
+
+    @tex.setter
+    def tex(self, t):
+        self._tex = t if t is not None else r'{\rm %s}' % self.name
+
+    @property
+    def tag(self):
+        return self._tag
+
+    @tag.setter
+    def tag(self, t):
+        if t is None: self._tag = ParamTag.NONE
+        else:
+            assert t in ParamTag
+            self._tag = t
+
+
+class ParamSet(Sequence):
+    """Container class for a set of parameters."""
+    def __init__(self, *args):
+        param_sequence = []
+        for arg in args:
+            try:
+                param_sequence.extend(arg)
+            except TypeError:
+                param_sequence.append(arg)
+
+        if len(param_sequence) != 0:
+            # Disallow duplicated params
+            all_names = [p.name for p in param_sequence]
+            unique_names = set(all_names)
+            if len(unique_names) != len(all_names):
+                duplicates = set([x for x in all_names if all_names.count(x) > 1])
+                raise ValueError('Duplicate definitions found for param(s): ' +
+                                 ', '.join(str(e) for e in duplicates))
+
+        # Elements of list must be Param type
+        assert all([isinstance(x, Param) for x in param_sequence]), \
+                'All params must be of type "Param"'
+
+        self._params = param_sequence
+
+    def __len__(self):
+        return len(self._params)
+
+    def __getitem__(self, i):
+        if isinstance(i, int):
+            return self._params[i]
+        elif isinstance(i, basestring):
+            return self._by_name[i]
+
+    def __getattr__(self, attr):
+        try:
+            return super(ParamSet, self).__getattribute__(attr)
+        except AttributeError:
+            t, v, tb = sys.exc_info()
+            try:
+                return self[attr]
+            except KeyError:
+                raise t, v, tb
+
+    def __iter__(self):
+        return iter(self._params)
+
+    def __str__(self):
+        o = '\n'
+        for obj in self._params:
+            o += '== {0:<15} = {1:<15}, tag={2:<15}\n'.format(
+                obj.name, obj.value, obj.tag
+            )
+        return o
+
+    @property
+    def _by_name(self):
+        return {obj.name: obj for obj in self._params}
+
+    @property
+    def names(self):
+        return tuple([obj.name for obj in self._params])
+
+    @property
+    def labels(self):
+        return tuple([obj.tex for obj in self._params])
+
+    @property
+    def values(self):
+        return tuple([obj.value for obj in self._params])
+
+    @property
+    def seeds(self):
+        return tuple([obj.seed for obj in self._params])
+
+    @property
+    def ranges(self):
+        return tuple([obj.ranges for obj in self._params])
+
+    @property
+    def stds(self):
+        return tuple([obj.std for obj in self._params])
+
+    @property
+    def tags(self):
+        return tuple([obj.tag for obj in self._params])
+
+    @property
+    def params(self):
+        return self._params
+
+    def to_dict(self):
+        return {obj.name: obj.value for obj in self._params}
+
+    def from_tag(self, tag, values=False, index=False, invert=False):
+        if values and index: assert 0
+        tag = np.atleast_1d(tag)
+        if not invert:
+            ps = [(idx, obj) for idx, obj in enumerate(self._params)
+                  if obj.tag in tag]
+        else:
+            ps = [(idx, obj) for idx, obj in enumerate(self._params)
+                  if obj.tag not in tag]
+        if values:
+            return tuple([io[1].value for io in ps])
+        elif index:
+            return tuple([io[0] for io in ps])
+        else:
+            return ParamSet([io[1] for io in ps])
+
+
+def get_paramsets(args):
+    """Make the paramsets for generating the Asmimov MC sample and also running
+    the MCMC.
+    """
+    asimov_paramset = []
+    llh_paramset = []
+    if args.likelihood == Likelihood.GOLEMFIT:
+        nuisance_paramset = define_nuisance()
+        asimov_paramset.extend(nuisance_paramset.params)
+        llh_paramset.extend(nuisance_paramset.params)
+        for parm in nuisance_paramset:
+            parm.value = args.__getattribute__(parm.name)
+    tag = ParamTag.BESTFIT
+    flavour_angles = fr_to_angles(args.measured_ratio)
+    asimov_paramset.extend([
+        Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag),
+        Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag),
+    ])
+    asimov_paramset = ParamSet(asimov_paramset)
+
+    if not args.fix_mixing and not args.fix_mixing_almost:
+        tag = ParamTag.MMANGLES
+        e = 1e-10
+        llh_paramset.extend([
+            Param(name='s_12^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag),
+            Param(name='c_13^4', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag),
+            Param(name='s_23^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag),
+            Param(name='dcp', value=np.pi, ranges=[0.+e, 2*np.pi-e], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag)
+        ])
+    if args.fix_mixing_almost:
+        tag = ParamTag.MMANGLES
+        llh_paramset.extend([
+            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:
+        tag = ParamTag.SCALE
+        llh_paramset.append(
+            Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
+                  tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag)
+        )
+    if not args.fix_source_ratio:
+        tag = ParamTag.SRCANGLES
+        llh_paramset.extend([
+            Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag),
+            Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag)
+        ])
+    llh_paramset = ParamSet(llh_paramset)
+    return asimov_paramset, llh_paramset
diff --git a/utils/plot.py b/utils/plot.py
index 0b82675..66c888c 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -86,26 +86,6 @@ def interval(arr, percentile=68.):
     return sarr[curr_low], center, sarr[curr_up]
 
 
-def plot_argparse(parser):
-    """Arguments for plotting."""
-    parser.add_argument(
-        '--plot-angles', type=misc_utils.parse_bool, default='True',
-        help='Plot MCMC triangle in the angles space'
-    )
-    parser.add_argument(
-        '--plot-elements', type=misc_utils.parse_bool, default='False',
-        help='Plot MCMC triangle in the mixing elements space'
-    )
-    parser.add_argument(
-        '--plot-bayes', type=misc_utils.parse_bool, default='False',
-        help='Plot Bayes factor'
-    )
-    parser.add_argument(
-        '--plot-angles-limit', type=misc_utils.parse_bool, default='False',
-        help='Plot limit vs BSM angles'
-    )
-
-
 def flat_angles_to_u(x):
     """Convert from angles to mixing elements."""
     return abs(angles_to_u(x)).astype(np.float32).flatten().tolist()
@@ -138,32 +118,19 @@ def gen_figtext(args):
     mr1, mr2, mr3 = args.measured_ratio
     if args.fix_source_ratio:
         sr1, sr2, sr3 = args.source_ratio
-        if args.fix_scale:
-            t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \
-                    'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \
-                    '{5:.2f}]\nDimension = {6}\nScale = {7}'.format(
-                        sr1, sr2, sr3, mr1, mr2, mr3, args.dimension,
-                        int(args.energy), args.scale
-                    )
-        else:
-            t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \
-                    'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \
-                    '{5:.2f}]\nDimension = {6}'.format(
-                        sr1, sr2, sr3, mr1, mr2, mr3, args.dimension,
-                        int(args.energy)
-                    )
+        t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \
+                'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \
+                '{5:.2f}]\nDimension = {6}'.format(
+                    sr1, sr2, sr3, mr1, mr2, mr3, args.dimension,
+                    int(args.energy)
+                )
     else:
-        if args.fix_scale:
-            t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \
-                    '{2:.2f}]\nDimension = {3}\nScale = {4}'.format(
-                        mr1, mr2, mr3, args.dimension, int(args.energy),
-                        args.scale
-                    )
-        else:
-            t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \
-                    '{2:.2f}]\nDimension = {3}'.format(
-                        mr1, mr2, mr3, args.dimension, int(args.energy)
-                    )
+        t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \
+                '{2:.2f}]\nDimension = {3}'.format(
+                    mr1, mr2, mr3, args.dimension, int(args.energy)
+                )
+    if args.fix_scale:
+        t += 'Scale = {0}'.format(args.scale)
     if args.likelihood is Likelihood.GAUSSIAN:
         t += '\nSigma = {0:.3f}'.format(args.sigma_ratio)
     if args.energy_dependance is EnergyDependance.SPECTRAL:
@@ -176,7 +143,7 @@ def gen_figtext(args):
     return t
 
 
-def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
+def chainer_plot(infile, outfile, outformat, args, llh_paramset):
     """Make the triangle plot."""
     if not args.plot_angles and not args.plot_elements:
         return
@@ -187,8 +154,8 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
     misc_utils.make_dir(outfile)
     fig_text = gen_figtext(args)
 
-    axes_labels = mcmc_paramset.labels
-    ranges = mcmc_paramset.ranges
+    axes_labels = llh_paramset.labels
+    ranges = llh_paramset.ranges
 
     if args.plot_angles:
         print "Making triangle plots"
@@ -208,7 +175,7 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
                         ), fontsize=10)
 
         # if not args.fix_mixing:
-        #     sc_index = mcmc_paramset.from_tag(ParamTag.SCALE, index=True)
+        #     sc_index = llh_paramset.from_tag(ParamTag.SCALE, index=True)
         #     itv = interval(Tchain[:,sc_index], percentile=90.)
         #     mpl.pyplot.figtext(
         #         0.5, 0.3, 'Scale 90% Interval = [1E{0}, 1E{1}], Center = '
@@ -222,11 +189,11 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
         print "Making triangle plots"
         if args.fix_mixing_almost:
             raise NotImplementedError
-        nu_index = mcmc_paramset.from_tag(ParamTag.NUISANCE, index=True)
-        fr_index = mcmc_paramset.from_tag(ParamTag.MMANGLES, index=True)
-        sc_index = mcmc_paramset.from_tag(ParamTag.SCALE, index=True)
+        nu_index = llh_paramset.from_tag(ParamTag.NUISANCE, index=True)
+        fr_index = llh_paramset.from_tag(ParamTag.MMANGLES, index=True)
+        sc_index = llh_paramset.from_tag(ParamTag.SCALE, index=True)
         if not args.fix_source_ratio:
-            sr_index = mcmc_paramset.from_tag(ParamTag.SRCANGLES, index=True)
+            sr_index = llh_paramset.from_tag(ParamTag.SRCANGLES, index=True)
 
         nu_elements = raw[:,nu_index]
         fr_elements = np.array(map(flat_angles_to_u, raw[:,fr_index]))
@@ -264,45 +231,39 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
             g.export(outfile+'_elements.'+of)
 
 
-def bayes_factor_plot(dirname, outfile, outformat, args):
-    """Make Bayes factor plot."""
-    if not args.plot_bayes: return
-    print "Making Bayes Factor plot"
-    print 'dirname', dirname
+def plot_multinest(data, outfile, outformat, args, scale_param, label=None):
+    """Make MultiNest factor plot."""
+    print "Making MultiNest Factor plot"
     fig_text = gen_figtext(args)
+    if label is not None: fig_text += '\n' + label
 
-    raw = []
-    for root, dirs, filenames in os.walk(dirname):
-        for fn in filenames:
-            if fn[-4:] == '.npy':
-                raw.append(np.load(os.path.join(root, fn)))
-    raw = np.sort(np.vstack(raw), axis=0)
-    print 'raw', raw
-    print 'raw.shape', raw.shape
-    scales, evidences = raw.T
-    null = evidences[0]
+    print 'data.shape', data.shape
+    scales, evidences = data.T
+    min_idx = np.argmin(scales)
+    null = evidences[min_idx]
 
     reduced_ev = -(evidences - null)
 
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
 
-    ax.set_xlim(np.log10(args.scale_region))
-    ax.set_xlabel(r'${\rm log}_{10} \Lambda ' + get_units(args.dimension) +r'$')
+    ax.set_xlim(scale_param.ranges)
+    ax.set_xlabel('$'+scale_param.tex+'$')
     ax.set_ylabel(r'Bayes Factor')
 
     ax.plot(scales, reduced_ev)
 
     for ymaj in ax.yaxis.get_majorticklocs():
-        ax.axhline(y=ymaj, ls='-', color='gray', alpha=0.4, linewidth=1)
+        ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.3, linewidth=1)
     for xmaj in ax.xaxis.get_majorticklocs():
-        ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1)
+        ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.3, linewidth=1)
 
     at = AnchoredText(
-        fig_text, prop=dict(size=7), frameon=True, loc=2
+        '\n'+fig_text, prop=dict(size=7), frameon=True, loc=2
     )
-    at.patch.set_boxstyle("round,pad=0.,rounding_size=0.5")
+    at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5")
     ax.add_artist(at)
+    misc_utils.make_dir(outfile)
     for of in outformat:
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
 
@@ -326,7 +287,8 @@ def plot_BSM_angles_limit(dirname, outfile, outformat, args, bayesian):
         for fn in filenames:
             if fn[-4:] == '.npy':
                 raw.append(np.load(os.path.join(root, fn)))
-    raw = np.sort(np.vstack(raw), axis=0)
+    raw = np.vstack(raw)
+    raw = raw[np.argsort(raw[:,0])]
     print 'raw', raw
     print 'raw.shape', raw.shape
     sc_ranges = (