1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
|
# 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 import llh as llh_utils
from utils.misc import gen_identifier, make_dir, parse_bool, solve_ratio
def CubePrior(cube, ndim, n_params):
pass
def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
args):
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.ranges
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
llh = llh_utils.ln_prob(
theta=theta,
args=args,
asimov_paramset=asimov_paramset,
llh_paramset=llh_paramset,
)
return llh
def mn_argparse(parser):
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(
'--run-mn', type=parse_bool, default='True',
help='Run MultiNest'
)
def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, prefix='mn'):
"""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,
)
if args.run_mn:
make_dir(prefix)
print 'Running evidence calculation for {0}'.format(prefix)
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,
resume = True,
verbose = True
)
analyser = analyse.Analyzer(
outputfiles_basename=prefix, n_params=n_params
)
return analyser.get_stats()['global evidence']
|
