|
| 1 | +# Copyright (c) 2015-2020 by the parties listed in the AUTHORS file. |
| 2 | +# All rights reserved. Use of this source code is governed by |
| 3 | +# a BSD-style license that can be found in the LICENSE file. |
| 4 | + |
| 5 | +import numpy as np |
| 6 | +import traitlets |
| 7 | +from astropy import units as u |
| 8 | + |
| 9 | +from ..accelerator import ImplementationType |
| 10 | +from ..observation import default_values as defaults |
| 11 | +from ..timing import function_timer |
| 12 | +from ..traits import Bool, Instance, Int, Unicode, trait_docs |
| 13 | +from ..utils import Environment, Logger |
| 14 | +from ..qarray import mult, to_iso_angles |
| 15 | +from .operator import Operator |
| 16 | + |
| 17 | +@trait_docs |
| 18 | +class DerivativesWeights(Operator): |
| 19 | + """Operator which generates pointing weights for I and derivatives of I with |
| 20 | + respect to theta and phi, to order 1 (mode dI) or 2 (mode d2I). |
| 21 | +
|
| 22 | + Given the individual detector pointing, this computes the pointing weights |
| 23 | + assuming that the detector is a linear polarizer followed by a total |
| 24 | + power measurement. By definition, the detector coordinate frame has the X-axis |
| 25 | + aligned with the polarization sensitive direction. An optional dictionary of |
| 26 | + beam error factors may be specified for each observation. |
| 27 | + |
| 28 | + These factors are an overall calibration factor cal, beam centroid error dx/dy, |
| 29 | + differential beam fwhm dsigma, and ellipticity dp/dc. Since we are focused |
| 30 | + on total intensity, there is no HWP term or detector polarisation efficiency. |
| 31 | +
|
| 32 | + The timestream model without a HWP in COSMO convention is: |
| 33 | +
|
| 34 | + .. math:: |
| 35 | + d = cal*I + d_\\theta I \\left[ dx\\sin\\psi - dy\\cos\\psi \\right] |
| 36 | + + d_\\phi I \\left[ -dx\\cos\\psi - dy\\sin\\psi + (dp\\sin(2\\psi) - dc\\cos(2\\psi))\frac{\\cos\\theta}{\\sin\\theta} \\right] |
| 37 | + + d^2_\\theta I \\left[dsigma + dp\\cos(2\\psi) - dc\\sin(2\\psi)\\right] |
| 38 | + + d_\\phi d_\\theta I \\left[-2dp\\sin(2\\psi) + 2dc\\cos(2\\psi) \\right] |
| 39 | + + d^2_\\phi I \\left[\\ dsigma + dp\\cos(2\\psi) + dc\\sin(2\\psi) \\right] |
| 40 | +
|
| 41 | + The detector orientation angle "psi" in COSMO convention is measured in a |
| 42 | + right-handed sense from the local meridian. |
| 43 | +
|
| 44 | + By default, this operator uses the "COSMO" convention for Q/U. If the "IAU" trait |
| 45 | + is set to True, then resulting weights will differ as psi will jump around. |
| 46 | +
|
| 47 | + If the view trait is not specified, then this operator will use the same data |
| 48 | + view as the detector pointing operator when computing the pointing matrix pixels |
| 49 | + and weights. |
| 50 | +
|
| 51 | + """ |
| 52 | + |
| 53 | + # Class traits |
| 54 | + |
| 55 | + API = Int(0, help="Internal interface version for this operator") |
| 56 | + |
| 57 | + detector_pointing = Instance( |
| 58 | + klass=Operator, |
| 59 | + allow_none=True, |
| 60 | + help="Operator that translates boresight pointing into detector frame", |
| 61 | + ) |
| 62 | + |
| 63 | + mode = Unicode("dI", help="The Stokes weights to generate (dI or d2I)") |
| 64 | + |
| 65 | + view = Unicode( |
| 66 | + None, allow_none=True, help="Use this view of the data in all observations" |
| 67 | + ) |
| 68 | + |
| 69 | + weights = Unicode( |
| 70 | + defaults.weights, help="Observation detdata key for output weights" |
| 71 | + ) |
| 72 | + |
| 73 | + single_precision = Bool(False, help="If True, use 32bit float in output") |
| 74 | + |
| 75 | + IAU = Bool(False, help="If True, use the IAU convention rather than COSMO") |
| 76 | + |
| 77 | + @traitlets.validate("detector_pointing") |
| 78 | + def _check_detector_pointing(self, proposal): |
| 79 | + detpointing = proposal["value"] |
| 80 | + if detpointing is not None: |
| 81 | + if not isinstance(detpointing, Operator): |
| 82 | + raise traitlets.TraitError( |
| 83 | + "detector_pointing should be an Operator instance" |
| 84 | + ) |
| 85 | + # Check that this operator has the traits we expect |
| 86 | + for trt in [ |
| 87 | + "view", |
| 88 | + "boresight", |
| 89 | + "shared_flags", |
| 90 | + "shared_flag_mask", |
| 91 | + "det_mask", |
| 92 | + "quats", |
| 93 | + "coord_in", |
| 94 | + "coord_out", |
| 95 | + ]: |
| 96 | + if not detpointing.has_trait(trt): |
| 97 | + msg = f"detector_pointing operator should have a '{trt}' trait" |
| 98 | + raise traitlets.TraitError(msg) |
| 99 | + return detpointing |
| 100 | + |
| 101 | + @traitlets.validate("mode") |
| 102 | + def _check_mode(self, proposal): |
| 103 | + check = proposal["value"] |
| 104 | + if check not in ["dI", "d2I"]: |
| 105 | + raise traitlets.TraitError("Invalid mode (must be 'dI' or 'd2I')") |
| 106 | + return check |
| 107 | + |
| 108 | + def __init__(self, **kwargs): |
| 109 | + super().__init__(**kwargs) |
| 110 | + |
| 111 | + @function_timer |
| 112 | + def _exec(self, data, detectors=None, use_accel=None, **kwargs): |
| 113 | + env = Environment.get() |
| 114 | + log = Logger.get() |
| 115 | + if self.mode == "d2I": |
| 116 | + self._nnz = 6 |
| 117 | + else: |
| 118 | + self._nnz = 3 |
| 119 | + |
| 120 | + # Kernel selection |
| 121 | + implementation, use_accel = self.select_kernels(use_accel=use_accel) |
| 122 | + |
| 123 | + if self.detector_pointing is None: |
| 124 | + raise RuntimeError("The detector_pointing trait must be set") |
| 125 | + |
| 126 | + # Expand detector pointing |
| 127 | + quats_name = self.detector_pointing.quats |
| 128 | + |
| 129 | + view = self.view |
| 130 | + if view is None: |
| 131 | + # Use the same data view as detector pointing |
| 132 | + view = self.detector_pointing.view |
| 133 | + |
| 134 | + # Expand detector pointing |
| 135 | + self.detector_pointing.apply(data, detectors=detectors, use_accel=use_accel) |
| 136 | + |
| 137 | + for ob in data.obs: |
| 138 | + # Get the detectors we are using for this observation |
| 139 | + dets = ob.select_local_detectors( |
| 140 | + detectors, flagmask=self.detector_pointing.det_mask |
| 141 | + ) |
| 142 | + if len(dets) == 0: |
| 143 | + # Nothing to do for this observation |
| 144 | + continue |
| 145 | + |
| 146 | + # Check that our view is fully covered by detector pointing. If the |
| 147 | + # detector_pointing view is None, then it has all samples. If our own |
| 148 | + # view was None, then it would have been set to the detector_pointing |
| 149 | + # view above. |
| 150 | + if (view is not None) and (self.detector_pointing.view is not None): |
| 151 | + if ob.intervals[view] != ob.intervals[self.detector_pointing.view]: |
| 152 | + # We need to check intersection |
| 153 | + intervals = ob.intervals[self.view] |
| 154 | + detector_intervals = ob.intervals[self.detector_pointing.view] |
| 155 | + intersection = detector_intervals & intervals |
| 156 | + if intersection != intervals: |
| 157 | + msg = ( |
| 158 | + f"view {self.view} is not fully covered by valid " |
| 159 | + "detector pointing" |
| 160 | + ) |
| 161 | + raise RuntimeError(msg) |
| 162 | + |
| 163 | + # Create (or re-use) output data for the weights |
| 164 | + if self.single_precision: |
| 165 | + exists = ob.detdata.ensure( |
| 166 | + self.weights, |
| 167 | + sample_shape=(self._nnz,), |
| 168 | + dtype=np.float32, |
| 169 | + detectors=dets, |
| 170 | + accel=use_accel, |
| 171 | + ) |
| 172 | + else: |
| 173 | + exists = ob.detdata.ensure( |
| 174 | + self.weights, |
| 175 | + sample_shape=(self._nnz,), |
| 176 | + dtype=np.float64, |
| 177 | + detectors=dets, |
| 178 | + accel=use_accel, |
| 179 | + ) |
| 180 | + |
| 181 | + # Do we already have pointing for all requested detectors? |
| 182 | + if exists: |
| 183 | + # Yes |
| 184 | + if data.comm.group_rank == 0: |
| 185 | + msg = ( |
| 186 | + f"Group {data.comm.group}, ob {ob.name}, derivative weights " |
| 187 | + f"already computed for {dets}" |
| 188 | + ) |
| 189 | + log.verbose(msg) |
| 190 | + continue |
| 191 | + |
| 192 | + # FIXME: temporary hack until instrument classes are also pre-staged |
| 193 | + # to GPU |
| 194 | + focalplane = ob.telescope.focalplane |
| 195 | + focalplane.detector_data |
| 196 | + #Get the boresight pointing |
| 197 | + qbore = ob.shared["boresight_radec"] |
| 198 | + nsamp = len(qbore) |
| 199 | + ndets = len(dets) |
| 200 | + theta = np.empty(nsamp) |
| 201 | + psi = np.empty(nsamp) |
| 202 | + # Get the per-detector pointing for orientation/sine theta purposes |
| 203 | + for idet, d in enumerate(dets): |
| 204 | + theta, _, psi = to_iso_angles(mult(qbore, focalplane[d]["quat"])) |
| 205 | + psi -= focalplane[d]["pol_angle"].value |
| 206 | + wc = np.cos(psi) |
| 207 | + wc2 = np.cos(2*psi) |
| 208 | + ws = np.sin(psi) |
| 209 | + ws2 = np.sin(2*psi) |
| 210 | + inv_tan_theta = np.cos(theta)/np.sin(theta) |
| 211 | + # Get the per-detector calibration. For now we sidestep with the Nones |
| 212 | + cal = focalplane[d].get("cal", 1.0) |
| 213 | + fwhm = focalplane[d].get("FWHM", 1.4*u.arcmin).to(u.rad).value |
| 214 | + dx = focalplane[d].get("dx", 0.0) |
| 215 | + dy = focalplane[d].get("dy", 0.0) |
| 216 | + dsigma = focalplane[d].get("dsigma", 0.0) |
| 217 | + dp = focalplane[d].get("dp", 0.0) |
| 218 | + dc = focalplane[d].get("dc", 0.0) |
| 219 | + |
| 220 | + b_std = fwhm/np.sqrt(8*np.log(2)) #beam standard deviation |
| 221 | + |
| 222 | + weights = np.empty((nsamp, self._nnz)) |
| 223 | + weights[:,0] = cal # gain error |
| 224 | + weights[:,1] = dx * ws - dy * wc #dtheta |
| 225 | + weights[:,2] = -dx * wc - dy * ws + b_std**2 * (dp * ws2 - dc * wc2) * inv_tan_theta #dphi |
| 226 | + if self.mode == "d2I": |
| 227 | + weights[:,3] = b_std * dsigma + 0.5 * b_std * b_std * (dp * wc2 - dc * ws2) #d2theta |
| 228 | + weights[:,4] = b_std**2 * (-2.0 * dp * ws2 + 2.0 * dc * wc2) #dphi dtheta |
| 229 | + weights[:,5] = b_std * dsigma + b_std**2 * (dp * wc2 + dc * ws2) #dphi2 |
| 230 | + ob.detdata[self.weights][d, :] = weights |
| 231 | + return |
| 232 | + """ |
| 233 | + # Get the per-detector calibration |
| 234 | + if self.cal is None: |
| 235 | + cal = np.array([1.0 for x in dets], np.float64) |
| 236 | + else: |
| 237 | + cal = np.array([ob[self.cal][x] for x in dets], np.float64) |
| 238 | + cal = np.stack([cal for _ in range(nsamp)], axis=1) |
| 239 | + # Per-detector pointing error |
| 240 | + if self.dx is None: |
| 241 | + dx = np.array([0.0 for x in dets], np.float64) |
| 242 | + else: |
| 243 | + dx = np.array([ob[self.dx][x] for x in dets], np.float64) |
| 244 | + if self.dy is None: |
| 245 | + dy = np.array([0.0 for x in dets], np.float64) |
| 246 | + else: |
| 247 | + dy = np.array([ob[self.dy][x] for x in dets], np.float64) |
| 248 | + dx = np.stack([dx for _ in range(nsamp)], axis=1) |
| 249 | + dy = np.stack([dy for _ in range(nsamp)], axis=1) |
| 250 | + #Per-detector fwhm/sigma error |
| 251 | + if self.dsigma is None: |
| 252 | + dsigma = np.array([0.0 for x in dets], np.float64) |
| 253 | + else: |
| 254 | + dsigma = np.array([ob[self.dsigma][x] for x in dets], np.float64) |
| 255 | + dsigma = np.stack([dsigma for _ in range(nsamp)], axis=1) |
| 256 | + #Per-detector ellipticity |
| 257 | + if self.dp is None: |
| 258 | + dp = np.array([0.0 for x in dets], np.float64) |
| 259 | + else: |
| 260 | + dp = np.array([ob[self.dp][x] for x in dets], np.float64) |
| 261 | + if self.dc is None: |
| 262 | + dc = np.array([0.0 for x in dets], np.float64) |
| 263 | + else: |
| 264 | + dc = np.array([ob[self.dc][x] for x in dets], np.float64) |
| 265 | + dp = np.stack([dp for _ in range(nsamp)], axis=1) |
| 266 | + dc = np.stack([dc for _ in range(nsamp)], axis=1) |
| 267 | +
|
| 268 | + wc = np.cos(psi) |
| 269 | + wc2 = np.cos(2*psi) |
| 270 | + ws = np.sin(psi) |
| 271 | + ws2 = np.sin(2*psi) |
| 272 | + inv_tan_theta = np.cos(theta)/np.sin(theta) |
| 273 | +
|
| 274 | + weights = np.empty((ndets,nsamp,self._nnz)) |
| 275 | + weights[:,:,0] = cal # gain error |
| 276 | + weights[:,:,1] = dx * ws - dy * wc #dtheta |
| 277 | + weights[:,:,2] = -dx * wc - dy * ws + (dp * ws2 - dc * wc2) * inv_tan_theta #dphi |
| 278 | + if self.mode == "d2I": |
| 279 | + weights[:,:,3] = dsigma + dp * wc2 - dc * ws2 #d2theta |
| 280 | + weights[:,:,4] = -2.0 * dp * ws2 + 2.0 * dc * wc2 #dphi dtheta |
| 281 | + weights[:,:,5] = dsigma + dp * wc2 + dc * ws2 #dphi2 |
| 282 | +
|
| 283 | + for idet, d in enumerate(dets): |
| 284 | + ob.detdata[self.weights][d, :] = weights[idet] |
| 285 | + """ |
| 286 | + |
| 287 | + |
| 288 | + def _finalize(self, data, **kwargs): |
| 289 | + return |
| 290 | + |
| 291 | + def _requires(self): |
| 292 | + req = self.detector_pointing.requires() |
| 293 | + if "detdata" not in req: |
| 294 | + req["detdata"] = list() |
| 295 | + req["detdata"].append(self.weights) |
| 296 | + if self.cal is not None: |
| 297 | + req["meta"].append(self.cal) |
| 298 | + if self.hwp_angle is not None: |
| 299 | + req["shared"].append(self.hwp_angle) |
| 300 | + if self.view is not None: |
| 301 | + req["intervals"].append(self.view) |
| 302 | + return req |
| 303 | + |
| 304 | + def _provides(self): |
| 305 | + prov = self.detector_pointing.provides() |
| 306 | + prov["detdata"].append(self.weights) |
| 307 | + return prov |
| 308 | + |
| 309 | + def _implementations(self): |
| 310 | + return [ |
| 311 | + ImplementationType.DEFAULT, |
| 312 | + ImplementationType.COMPILED, |
| 313 | + ImplementationType.NUMPY, |
| 314 | + ImplementationType.JAX, |
| 315 | + ] |
| 316 | + |
| 317 | + def _supports_accel(self): |
| 318 | + if (self.detector_pointing is not None) and ( |
| 319 | + self.detector_pointing.supports_accel() |
| 320 | + ): |
| 321 | + return True |
| 322 | + else: |
| 323 | + return False |
0 commit comments