Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInference 4.1.9.1-8a6b96f
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
cbcBayesThermoInt.py
Go to the documentation of this file.
1##python
2
3from optparse import OptionParser
4import sys
5import os
6import matplotlib
7matplotlib.use('Agg')
8import matplotlib.pyplot as pp
9import numpy as np
10try:
11 from scipy.integrate import trapezoid
12except ImportError:
13 # FIXME: Remove this once we require scipy >=1.6.0.
14 from scipy.integrate import trapz as trapezoid
15
16matplotlib.rcParams['text.usetex'] = True
17
18def extract_temp(filename):
19 """Extracts the PTMCMC temperature from the header lines of the
20 given file."""
21 with open(filename, 'r') as inp:
22 for line in inp:
23 line=line.lstrip().split()
24 try:
25 i=line.index('Tchain')
26 #####################################################
27 # WARNING: hardcoded off-by-one adjustment to account
28 # for 'null likelihood' header name that splits into
29 # two list elements
30 #####################################################
31 return float(inp.next().split()[i-1])
32 except ValueError:
33 pass
34
35 raise ValueError('extract_temp: did not find header line with \'Tchain\'')
36
37def get_mean_logl(filename):
38 """Returns the mean value of log(L) from the given filename,
39 excluding the first 50% of samples as burn-in."""
40 with open(filename, 'r') as inp:
41 skip=0
42 for line in inp:
43 line=line.lstrip().split()
44 skip+=1
45 try:
46 line.index('cycle')
47 break
48 except ValueError:
49 pass
50
51 col = line.index('logl')
52
53 data=np.loadtxt(filename, skiprows=skip, usecols=(col,))
54 N=data.shape[0]
55
56 return np.mean(data[N/2:])
57
58if __name__=='__main__':
59 # Custom usage and help message
60 usage = """%s [-h] [--plotfile FILE] [--evfile FILE] OUTPUT_FILE [OUTPUT_FILE ...]
61
62Thermodynamic integration on PTMCMC samples.
63
64positional arguments:
65 OUTPUT_FILE PTMCMC output files""" % (os.path.basename(sys.argv[0]))
66
67 parser = OptionParser(usage=usage)
68 parser.add_option('--plotfile', metavar='FILE', default='evidence-integrand.pdf', help='plot output file')
69 parser.add_option('--evfile', metavar='FILE', default='evidence.dat', help='evidence output file')
70
71 (options, args) = parser.parse_args()
72
73 # Make positional arguments required
74 if len(args)==0:
75 parser.error('Positional filename argument(s) are required')
76
77 betas = np.array([1.0/extract_temp(f) for f in args])
78 logls = np.array([get_mean_logl(f) for f in args])
79
80 inds = np.argsort(betas)[::-1]
81
82 betas = betas[inds]
83 logls = logls[inds]
84
85 # Now extend to infinite temperature by copying the last <log(L)>.
86 # This works as long as the chains have extended to high enough
87 # temperature to sample the prior.
88 ebetas = np.concatenate((betas, [0.0]))
89 elogls = np.concatenate((logls, [logls[-1]]))
90
91 ebetas2 = np.concatenate((betas[::2], [0.0]))
92 elogls2 = np.concatenate((logls[::2], [logls[::2][-1]]))
93
94 evidence = -trapezoid(elogls, ebetas)
95 evidence2 = -trapezoid(elogls2, ebetas2)
96
97 devidence = np.abs(evidence - evidence2)
98
99 pp.plot(betas, betas*logls)
100 pp.xscale('log')
101 pp.xlabel(r'$\beta$')
102 pp.ylabel(r'$\beta \left\langle \ln \mathcal{L} \right\rangle$')
103 pp.savefig(options.plotfile)
104
105 with open(options.evfile, 'w') as out:
106 out.write('# ln-Z delta-ln-Z\n')
107 out.write(str(evidence) + ' ' + str(devidence))
def get_mean_logl(filename)
Returns the mean value of log(L) from the given filename, excluding the first 50% of samples as burn-...
def extract_temp(filename)
Extracts the PTMCMC temperature from the header lines of the given file.