-
Notifications
You must be signed in to change notification settings - Fork 6
/
modifypdf.py
111 lines (86 loc) · 2.97 KB
/
modifypdf.py
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
102
103
104
105
106
107
108
109
110
111
"""
2014 February 12
Shane Bussmann
A collection of utilities for modifying posterior probability density
functions.
"""
from __future__ import print_function
from astropy.table import Table
import numpy
def trim(oldpdfloc, newpdfloc, niters=5000):
# get the last niters iterations
PDFdata = Table.read(oldpdfloc, format='ascii')
PDFdata = PDFdata[-niters:]
# write the trimmed list
PDFdata.write(newpdfloc, format='ascii')
def bigScoop(PDFdata, Nbuffer=300):
""" get model fits within Nbuffer of the best-fit model. """
lnprob = PDFdata['lnprob']
good = lnprob > lnprob.max() - Nbuffer
PDFdata = PDFdata[good]
# write the trimmed list
return PDFdata
def goodMu(PDFdata):
""" get model fits with mu < 100. """
# search for magnification measurements
lnprob = PDFdata['lnprob']
good = lnprob > lnprob.max() - 300
PDFdata = PDFdata[good]
# write the trimmed list
return PDFdata
def cleanColumns(PDFdata):
# get the last niters iterations
PDFkeys = PDFdata.keys()
for key in PDFkeys:
rms = numpy.std(PDFdata[key])
if rms == 0:
PDFdata.remove_column(key)
# write the trimmed list
return PDFdata
def prune(PDFdata, scaler=5.0, quiet=False):
# get the last niters iterations
#PDFdata = Table.read(oldpdfloc, format='ascii')
okok = PDFdata['lnprob'] * 0 == 0
PDFdata = PDFdata[okok]
if not quiet:
print("Pre-pruning, <Lnprob>: {:f}".format(PDFdata['lnprob'].mean()))
#import pdb; pdb.set_trace()
# identify the good fits
lnprob = PDFdata['lnprob']
#PDFdata['lnprob'] = lnprob
minlnprob = lnprob.max()
dlnprob = numpy.abs(lnprob - minlnprob)
medlnprob = numpy.median(dlnprob)
avglnprob = numpy.mean(dlnprob)
skewlnprob = numpy.abs(avglnprob - medlnprob)
rmslnprob = numpy.std(dlnprob)
inliers = (dlnprob < scaler*rmslnprob)
PDFdata = PDFdata[inliers]
medlnprob_previous = 0.
while skewlnprob > 0.1*medlnprob:
lnprob = PDFdata['lnprob']
minlnprob = lnprob.max()
dlnprob = numpy.abs(lnprob - minlnprob)
rmslnprob = numpy.std(dlnprob)
inliers = (dlnprob < scaler*rmslnprob)
PDFdatatmp = PDFdata[inliers]
if len(PDFdatatmp) == len(PDFdata):
inliers = (dlnprob < scaler/2.*rmslnprob)
PDFdata = PDFdata[inliers]
lnprob = PDFdata['lnprob']
dlnprob = numpy.abs(lnprob - minlnprob)
medlnprob = numpy.median(dlnprob)
avglnprob = numpy.mean(dlnprob)
skewlnprob = numpy.abs(avglnprob - medlnprob)
if not quiet:
print(medlnprob, avglnprob, skewlnprob)
if medlnprob == medlnprob_previous:
scaler /= 1.5
medlnprob_previous = medlnprob
#PDFdatagood = PDFdata.copy()
# identify the good fits
goodfits = (PDFdata['lnprob'] <= minlnprob)
PDFdata = PDFdata[goodfits]
# return the trimmed list
#PDFdata.write(newpdfloc, format='ascii')
return PDFdata