arbitrary precision calculations

8 files changed, 96 insertions, 52 deletions

diff --git a/plot_sens.py b/plot_sens.py
index de7b348..3d1181a 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -244,6 +244,7 @@ def main():
             if args.likelihood is Likelihood.GAUSSIAN:
                 infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
             infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/fr_stat'.format(
+            # infile += '/DIM{0}/fix_ifr/100TeV/{1}/{2}/{3}/fr_stat'.format(
                 dim, *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
             ) + misc_utils.gen_identifier(argsc)
             print '== {0:<25} = {1}'.format('infile', infile)
@@ -290,6 +291,7 @@ def main():
             if args.likelihood is Likelihood.GAUSSIAN:
                 base_infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
             base_infile += '/DIM{0}/fix_ifr'.format(dim)
+            # base_infile += '/DIM{0}/fix_ifr/100TeV'.format(dim)
 
             for isrc, src in enumerate(args.source_ratios):
                 argsc.source_ratio = src
diff --git a/sens.py b/sens.py
index 581a3ca..58f8231 100755
--- a/sens.py
+++ b/sens.py
@@ -25,7 +25,6 @@ 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
@@ -236,12 +235,12 @@ def main():
             if args.run_method in fixed_angle_categ:
                 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/4., dtype=DTYPE)**2
+                    mmangles[idx_scen].value = np.sin(np.pi/4.)**2
                     """s_12^2 or s_23^2"""
                     mmangles[1].value = 1.
                     """c_13^4"""
                 elif idx_scen == 1:
-                    mmangles[idx_scen].value = np.cos(np.pi/4., dtype=DTYPE)**4
+                    mmangles[idx_scen].value = np.cos(np.pi/4.)**4
                     """c_13^4"""
 
             for idx_an, an in enumerate(scan_angles):
@@ -274,18 +273,17 @@ def main():
                     print '|||| SCALE = {0:.0E}'.format(np.power(10, sc))
                     scale.value = sc
                     if args.stat_method is StatCateg.BAYESIAN:
-                        try:
-                            stat = mn_utils.mn_evidence(
-                                mn_paramset     = sens_paramset,
-                                llh_paramset    = llh_paramset,
-                                asimov_paramset = asimov_paramset,
-                                args            = args,
-                                fitter          = fitter
-                            )
-                        except:
+                        stat = mn_utils.mn_evidence(
+                            mn_paramset     = sens_paramset,
+                            llh_paramset    = llh_paramset,
+                            asimov_paramset = asimov_paramset,
+                            args            = args,
+                            fitter          = fitter
+                        )
+                        if stat is None:
                             print 'Failed run, continuing'
-                            # raise
-                            continue
+                            raise
+                            # continue
                         print '## Evidence = {0}'.format(stat)
                     elif args.stat_method is StatCateg.FREQUENTIST:
                         def fn(x):
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index 0652705..419ce2e 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -27,9 +27,9 @@ GLOBAL_PARAMS = {}
 sens_eval_bin = 'true' # set to 'all' to run normally
 GLOBAL_PARAMS.update(dict(
     sens_run      = 'True',
-    run_method    = 'corr_angle', # full, fixed_angle, corr_angle
+    run_method    = 'fixed_angle', # full, fixed_angle, corr_angle
     stat_method   = 'bayesian',
-    sens_bins     = 15,
+    sens_bins     = 20,
     seed          = 'None'
 ))
 
@@ -41,8 +41,8 @@ GLOBAL_PARAMS.update(dict(
 ))
 
 # FR
-dimension         = [3]
-# dimension         = [3, 6]
+# dimension         = [3]
+dimension         = [3, 6]
 GLOBAL_PARAMS.update(dict(
     threads           = 1,
     binning           = '6e4 1e7 20',
@@ -65,7 +65,7 @@ GLOBAL_PARAMS.update(dict(
 # GolemFit
 GLOBAL_PARAMS.update(dict(
     ast  = 'p2_0',
-    data = 'real'
+    data = 'asimov'
 ))
 
 # Plot
diff --git a/utils/fr.py b/utils/fr.py
index f98c730..51864a3 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -14,14 +14,61 @@ from functools import partial
 import numpy as np
 
 from utils.enums import EnergyDependance
-from utils.misc import DTYPE, CDTYPE, PI
 from utils.misc import enum_parse, parse_bool
 
+import mpmath as mp
+mp.mp.dps = 100 # Computation precision
+
+# DTYPE  = np.float128
+# CDTYPE = np.complex256
+# PI     = np.arccos(DTYPE(-1))
+# SQRT   = np.sqrt
+# COS    = np.cos
+# SIN    = np.sin
+# ACOS   = np.arccos
+# ASIN   = np.arcsin
+# EXP    = np.exp
+
+DTYPE  = mp.mpf
+CDTYPE = mp.mpc
+PI     = mp.pi
+SQRT   = mp.sqrt
+COS    = mp.cos
+SIN    = mp.sin
+ACOS   = mp.acos
+ASIN   = mp.asin
+EXP    = mp.exp
 
 MASS_EIGENVALUES = [7.40E-23, 2.515E-21]
 """SM mass eigenvalues"""
 
 
+def determinant(x):
+    """Calculate the determininant of a 3x3 matrix.
+
+    Parameters
+    ----------
+    x : ndarray, shape = (3, 3)
+
+    Returns
+    ----------
+    float determinant
+
+    Examples
+    ----------
+    >>> print determinant(
+    >>>     [[-1.65238188-0.59549718j,  0.27486548-0.18437467j, -1.35524534-0.38542072j],
+    >>>      [-1.07480906+0.29630449j, -0.47808456-0.80316821j, -0.88609356-1.50737308j],
+    >>>      [-0.14924144-0.99230446j,  0.49504234+0.63639805j, 2.29258915-0.36537507j]]
+    >>> )
+    (2.7797571563274688+3.0841795325804848j)
+
+    """
+    return (x[0][0] * (x[1][1] * x[2][2] - x[2][1] * x[1][2])
+           -x[1][0] * (x[0][1] * x[2][2] - x[2][1] * x[0][2])
+           +x[2][0] * (x[0][1] * x[1][2] - x[1][1] * x[0][2]))
+
+
 def angles_to_fr(src_angles):
     """Convert angular projection of the source flavour ratio back into the
     flavour ratio.
@@ -41,13 +88,13 @@ def angles_to_fr(src_angles):
     (0.38340579025361626, 0.16431676725154978, 0.45227744249483393)
 
     """
-    sphi4, c2psi = src_angles
+    sphi4, c2psi = map(DTYPE, src_angles)
 
-    psi = (0.5)*np.arccos(c2psi, dtype=DTYPE)
+    psi = (0.5)*ACOS(c2psi)
 
-    sphi2 = np.sqrt(sphi4, dtype=DTYPE)
+    sphi2 = SQRT(sphi4)
     cphi2 = 1. - sphi2
-    spsi2 = np.sin(psi, dtype=DTYPE)**2
+    spsi2 = SIN(psi)**2
     cspi2 = 1. - spsi2
 
     x = sphi2*cspi2
@@ -78,27 +125,27 @@ def angles_to_u(bsm_angles):
            [ 0.28614067-0.42427084j, -0.64749908-0.21213542j,  0.52331757+0.j        ]])
 
     """
-    s12_2, c13_4, s23_2, dcp = bsm_angles
+    s12_2, c13_4, s23_2, dcp = map(DTYPE, bsm_angles)
     dcp = CDTYPE(dcp)
 
     c12_2 = 1. - s12_2
-    c13_2 = np.sqrt(c13_4, dtype=DTYPE)
+    c13_2 = SQRT(c13_4)
     s13_2 = 1. - c13_2
     c23_2 = 1. - 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))
+    t12 = ASIN(SQRT(s12_2))
+    t13 = ACOS(SQRT(c13_2))
+    t23 = ASIN(SQRT(s23_2))
 
-    c12 = np.cos(t12)
-    s12 = np.sin(t12)
-    c13 = np.cos(t13)
-    s13 = np.sin(t13)
-    c23 = np.cos(t23)
-    s23 = np.sin(t23)
+    c12 = COS(t12)
+    s12 = SIN(t12)
+    c13 = COS(t13)
+    s13 = SIN(t13)
+    c23 = COS(t23)
+    s23 = SIN(t23)
 
     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)
+    p2 = np.array([[c13 , 0   , s13*EXP(-1j*dcp)] , [0    , 1   , 0]   , [-s13*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)
@@ -141,15 +188,15 @@ def cardano_eqn(ham):
 
     a = -np.trace(ham)
     b = DTYPE(1)/2 * ((np.trace(ham))**DTYPE(2) - np.trace(np.dot(ham, ham)))
-    c = -np.linalg.det(ham)
+    c = -determinant(ham)
 
     Q = (DTYPE(1)/9) * (a**DTYPE(2) - DTYPE(3)*b)
     R = (DTYPE(1)/54) * (DTYPE(2)*a**DTYPE(3) - DTYPE(9)*a*b + DTYPE(27)*c)
-    theta = np.arccos(R / np.sqrt(Q**DTYPE(3)))
+    theta = ACOS(R / SQRT(Q**DTYPE(3)))
 
-    E1 = -DTYPE(2) * np.sqrt(Q) * np.cos(theta/DTYPE(3)) - (DTYPE(1)/3)*a
-    E2 = -DTYPE(2) * np.sqrt(Q) * np.cos((theta - DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
-    E3 = -DTYPE(2) * np.sqrt(Q) * np.cos((theta + DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
+    E1 = -DTYPE(2) * SQRT(Q) * COS(theta/DTYPE(3)) - (DTYPE(1)/3)*a
+    E2 = -DTYPE(2) * SQRT(Q) * COS((theta - DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
+    E3 = -DTYPE(2) * SQRT(Q) * COS((theta + DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(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]
@@ -163,9 +210,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(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)
+    N1 = SQRT(np.abs(A1*B1)**2 + np.abs(A1*C1)**2 + np.abs(B1*C1)**2)
+    N2 = SQRT(np.abs(A2*B2)**2 + np.abs(A2*C2)**2 + np.abs(B2*C2)**2)
+    N3 = SQRT(np.abs(A3*B3)**2 + np.abs(A3*C3)**2 + np.abs(B3*C3)**2)
 
     mm = np.array([
         [np.conjugate(B1)*C1 / N1, np.conjugate(B2)*C2 / N2, np.conjugate(B3)*C3 / N3],
@@ -383,7 +430,7 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
     mass_matrix = np.array(
         [[0, 0, 0], [0, mass_eigenvalues[0], 0], [0, 0, mass_eigenvalues[1]]]
     )
-    sm_ham = (1./(2*energy))*np.matmul(sm_u, np.matmul(mass_matrix, sm_u.conj().T))
+    sm_ham = (1./(2*energy))*np.dot(sm_u, np.dot(mass_matrix, sm_u.conj().T))
     if no_bsm:
         eg_vector = cardano_eqn(sm_ham)
     else:
@@ -391,7 +438,7 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
         scale_matrix = np.array(
             [[0, 0, 0], [0, sc1, 0], [0, 0, sc2]]
         )
-        bsm_term = (energy**(dim-3)) * np.matmul(new_physics_u, np.matmul(scale_matrix, new_physics_u.conj().T))
+        bsm_term = (energy**(dim-3)) * np.dot(new_physics_u, np.dot(scale_matrix, new_physics_u.conj().T))
         bsm_ham = sm_ham + bsm_term
         eg_vector = cardano_eqn(bsm_ham)
 
diff --git a/utils/gf.py b/utils/gf.py
index c316ed2..79f1dd0 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -66,6 +66,7 @@ def steering_params(args):
     # params.years = [999]
 
     params.minFitEnergy = args.binning[0] # GeV
+    params.load_data_from_text_file = False
 
     return params
 
diff --git a/utils/misc.py b/utils/misc.py
index 7e40d13..837c145 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -22,11 +22,6 @@ 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 f3595f9..387ff90 100644
--- a/utils/multinest.py
+++ b/utils/multinest.py
@@ -95,6 +95,7 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
         resume                     = False,
         verbose                    = True
     )
+    if result is None: return None
 
     analyser = analyse.Analyzer(
         outputfiles_basename=prefix, n_params=n_params
diff --git a/utils/plot.py b/utils/plot.py
index 65635bc..b63126b 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -405,7 +405,7 @@ def plot_sens_fixed_angle(data, outfile, outformat, args):
                 print 'limit = {0}'.format(lim)
                 label = '{0} : {1} : {2}'.format(*misc_utils.solve_ratio(src))
                 if lim < yranges[0]: yranges[0] = lim
-                if lim > yranges[1]: yranges[1] = lim+arr_len+1
+                if lim > yranges[1]: yranges[1] = lim+arr_len+2
                 # if lim > yranges[1]: yranges[1] = lim
                 xoff = 0.15
                 line = plt.Line2D(