| 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.