1
|
|
|
import numpy as np |
|
|
|
|
2
|
|
|
import scipy.interpolate |
3
|
|
|
import scipy.optimize |
4
|
|
|
import pandas as pd |
5
|
|
|
|
6
|
|
|
_INTEGRAL_OUTER_RADIUS = 15.0 |
7
|
|
|
|
8
|
|
|
|
9
|
|
|
class InvalidPSFError(ValueError): |
|
|
|
|
10
|
|
|
|
11
|
|
|
pass |
12
|
|
|
|
13
|
|
|
|
14
|
|
|
class PSFWrapper(object): |
|
|
|
|
15
|
|
|
|
16
|
|
|
def __init__(self, xs, ys, brightness_interp_x=None, brightness_interp_y=None): |
17
|
|
|
|
18
|
|
|
self._xs = xs |
19
|
|
|
self._ys = ys |
20
|
|
|
|
21
|
|
|
self._psf_interpolated = scipy.interpolate.InterpolatedUnivariateSpline(xs, ys, k=2, |
22
|
|
|
ext='raise', |
23
|
|
|
check_finite=True) |
24
|
|
|
|
25
|
|
|
# Memorize the total integral (will use it for normalization) |
26
|
|
|
|
27
|
|
|
self._total_integral = self._psf_interpolated.integral(self._xs[0], _INTEGRAL_OUTER_RADIUS) |
28
|
|
|
|
29
|
|
|
# Now compute the truncation radius, which is a very conservative measurement |
30
|
|
|
# of the size of the PSF |
31
|
|
|
|
32
|
|
|
self._truncation_radius = self.find_eef_radius(0.9999) |
33
|
|
|
|
34
|
|
|
# Let's also compute another measurement more appropriate for convolution |
35
|
|
|
self._kernel_radius = self.find_eef_radius(0.999) |
36
|
|
|
|
37
|
|
|
assert self._kernel_radius <= self._truncation_radius |
38
|
|
|
assert self._truncation_radius <= _INTEGRAL_OUTER_RADIUS |
39
|
|
|
|
40
|
|
|
# Prepare brightness interpolation |
41
|
|
|
|
42
|
|
|
if brightness_interp_x is None: |
43
|
|
|
|
44
|
|
|
brightness_interp_x, brightness_interp_y = self._prepare_brightness_interpolation_points() |
45
|
|
|
|
46
|
|
|
self._brightness_interp_x = brightness_interp_x |
47
|
|
|
self._brightness_interp_y = brightness_interp_y |
48
|
|
|
|
49
|
|
|
self._brightness_interpolation = scipy.interpolate.InterpolatedUnivariateSpline(brightness_interp_x, |
50
|
|
|
brightness_interp_y, |
51
|
|
|
k=2, |
52
|
|
|
ext='extrapolate', |
53
|
|
|
check_finite=True) |
54
|
|
|
|
55
|
|
|
def _prepare_brightness_interpolation_points(self): |
56
|
|
|
|
57
|
|
|
# Get the centers of the bins |
58
|
|
|
interp_x = (self._xs[1:] + self._xs[:-1]) / 2.0 |
59
|
|
|
|
60
|
|
|
# Compute the density |
61
|
|
|
|
62
|
|
|
interp_y = np.array(map(lambda (a, b): self.integral(a, b) / (np.pi * (b ** 2 - a ** 2)) / self._total_integral, |
|
|
|
|
63
|
|
|
zip(self._xs[:-1], self._xs[1:]))) |
64
|
|
|
|
65
|
|
|
# Add zero at r = _INTEGRAL_OUTER_RADIUS so that the extrapolated values will be correct |
66
|
|
|
interp_x = np.append(interp_x, [_INTEGRAL_OUTER_RADIUS]) |
67
|
|
|
interp_y = np.append(interp_y, [0.0]) |
68
|
|
|
|
69
|
|
|
return interp_x, interp_y |
70
|
|
|
|
71
|
|
|
def find_eef_radius(self, fraction): |
|
|
|
|
72
|
|
|
|
73
|
|
|
f = lambda r: fraction - self.integral(1e-4, r) / self._total_integral |
|
|
|
|
74
|
|
|
|
75
|
|
|
radius, status = scipy.optimize.brentq(f, 0.005, _INTEGRAL_OUTER_RADIUS, full_output = True) |
|
|
|
|
76
|
|
|
|
77
|
|
|
assert status.converged, "Brentq did not converged" |
78
|
|
|
|
79
|
|
|
return radius |
80
|
|
|
|
81
|
|
|
def brightness(self, r): |
|
|
|
|
82
|
|
|
|
83
|
|
|
return self._brightness_interpolation(r) |
84
|
|
|
|
85
|
|
|
@property |
86
|
|
|
def xs(self): |
|
|
|
|
87
|
|
|
""" |
88
|
|
|
X of the interpolation data |
89
|
|
|
""" |
90
|
|
|
return self._xs |
91
|
|
|
|
92
|
|
|
@property |
93
|
|
|
def ys(self): |
|
|
|
|
94
|
|
|
""" |
95
|
|
|
Y of the interpolation data |
96
|
|
|
""" |
97
|
|
|
return self._ys |
98
|
|
|
|
99
|
|
|
def combine_with_other_psf(self, other_psf, w1, w2): |
|
|
|
|
100
|
|
|
""" |
101
|
|
|
Return a PSF which is the linear interpolation between this one and the other one provided |
102
|
|
|
|
103
|
|
|
:param other_psf: another psf |
104
|
|
|
:param w1: weight for self (i.e., this PSF) |
105
|
|
|
:param w2: weight for the other psf |
106
|
|
|
:return: another PSF instance |
107
|
|
|
""" |
108
|
|
|
|
109
|
|
|
if isinstance(self, InvalidPSF) or isinstance(other_psf, InvalidPSF): |
110
|
|
|
return InvalidPSF() |
111
|
|
|
|
112
|
|
|
# Weight the ys |
113
|
|
|
new_ys = w1 * self.ys + w2 * other_psf.ys |
114
|
|
|
|
115
|
|
|
# Also weight the brightness interpolation points |
116
|
|
|
new_br_interp_y = w1 * self._brightness_interp_y + w2 * other_psf._brightness_interp_y |
|
|
|
|
117
|
|
|
|
118
|
|
|
return PSFWrapper(self.xs, new_ys, |
119
|
|
|
brightness_interp_x=self._brightness_interp_x, |
120
|
|
|
brightness_interp_y=new_br_interp_y) |
121
|
|
|
|
122
|
|
|
def to_pandas(self): |
|
|
|
|
123
|
|
|
|
124
|
|
|
items = (('xs', self._xs), ('ys', self._ys)) |
125
|
|
|
|
126
|
|
|
return pd.DataFrame.from_dict(dict(items)) |
127
|
|
|
|
128
|
|
|
@classmethod |
129
|
|
|
def from_pandas(cls, df): |
|
|
|
|
130
|
|
|
|
131
|
|
|
return cls(df.loc[:, 'xs'].values, df.loc[:, 'ys'].values) |
132
|
|
|
|
133
|
|
|
@classmethod |
134
|
|
|
def from_TF1(cls, tf1_instance): |
|
|
|
|
135
|
|
|
|
136
|
|
|
# Annoyingly, some PSFs for some Dec bins (at large Zenith angles) have |
137
|
|
|
# zero or negative integrals, i.e., they are useless. Return an unusable PSF |
138
|
|
|
# object in that case |
139
|
|
|
if tf1_instance.Integral(0, _INTEGRAL_OUTER_RADIUS) <= 0.0: |
140
|
|
|
|
141
|
|
|
return InvalidPSF() |
142
|
|
|
|
143
|
|
|
# Make interpolation |
144
|
|
|
xs = np.logspace(-3, np.log10(_INTEGRAL_OUTER_RADIUS), 500) |
|
|
|
|
145
|
|
|
ys = np.array(map(lambda x:tf1_instance.Eval(x), xs), float) |
|
|
|
|
146
|
|
|
|
147
|
|
|
assert np.all(np.isfinite(xs)) |
148
|
|
|
assert np.all(np.isfinite(ys)) |
149
|
|
|
|
150
|
|
|
instance = cls(xs, ys) |
151
|
|
|
|
152
|
|
|
instance._tf1 = tf1_instance.Clone() |
|
|
|
|
153
|
|
|
|
154
|
|
|
return instance |
155
|
|
|
|
156
|
|
|
def integral(self, a, b): |
|
|
|
|
157
|
|
|
|
158
|
|
|
return self._psf_interpolated.integral(a, b) |
159
|
|
|
|
160
|
|
|
@property |
161
|
|
|
def truncation_radius(self): |
|
|
|
|
162
|
|
|
return self._truncation_radius |
163
|
|
|
|
164
|
|
|
@property |
165
|
|
|
def total_integral(self): |
|
|
|
|
166
|
|
|
return self._total_integral |
167
|
|
|
|
168
|
|
|
@property |
169
|
|
|
def kernel_radius(self): |
|
|
|
|
170
|
|
|
return self._kernel_radius |
171
|
|
|
|
172
|
|
|
|
173
|
|
|
# This is a class that, whatever you try to use it for, will raise an exception. |
174
|
|
|
# This is to make sure that we never use an invalid PSF without knowing it |
175
|
|
|
class InvalidPSF(object): |
|
|
|
|
176
|
|
|
|
177
|
|
|
# It can be useful to copy an invalid PSF. For instance HAL.get_simulated_dataset() makes a |
178
|
|
|
# copy of the HAL instance, including detector response, which can contain InvalidPSF (which |
179
|
|
|
# is fine as long as they are not used). |
180
|
|
|
def __deepcopy__(self, memo): |
181
|
|
|
return InvalidPSF() |
182
|
|
|
|
183
|
|
|
def __getattribute__(self, item): |
184
|
|
|
|
185
|
|
|
if item == "__deepcopy__": |
186
|
|
|
return object.__getattribute__(self, item) |
187
|
|
|
raise InvalidPSFError("Trying to use an invalid PSF") |
188
|
|
|
|
The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:
If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.