From 360f073891225ffe913d86cffa253db22eb49ba8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Tue, 18 Jun 2024 05:21:20 -0400 Subject: [PATCH 01/13] Plans for arbitrary detector response stuff --- src/Projection.cxx | 52 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/src/Projection.cxx b/src/Projection.cxx index 0a9ae344..ae7951de 100644 --- a/src/Projection.cxx +++ b/src/Projection.cxx @@ -121,6 +121,58 @@ inline bool isNone(const bp::object &pyo) // shape [1][ndet][nrange]. More complicated schemes could be used, but this // seems to work well enough. +// Arbitrary detector response support +// ----------------------------------- +// Useful to support detectors with polarization response < 1. Also occasionally +// useful with detectors with 0 total intensity response. Currently the detector +// response is only encoded indirectly through the detector angle on the sky, part +// of the per-detector quaternion. Hard to cram the rest of the response into this, +// so must be handled with a new argument. +// +// spin_proj_factors: cos(2psi),sin(2psi) → detector on-sky TQU response +// psi is here the projected detector angle, which is computed by Pointer::GetCoords +// It makes sense for this to deal with angles since it's the result of a quaternion +// pointing calculation. +// +// Probably best to leave GetCoords alone, and instead augment spin_proj_factors with +// T and P support. For example: +// +// void spin_proj_factors(const double * coords, const double * response, FSIGNAL * projfacs) { +// const double c = coords[2]; +// const double s = coords[3]; +// projfacs[0] = response[0]; +// projfacs[1] = response[1]*(c*c - s*s); +// projfacs[2] = response[1]*(2*c*s); +// } +// +// None of the arguments to to_map and its relatives are a natural place to put response. +// It could be shoehorned into pofs by making it length 6, but the first four components +// would hold quaternion coefficients making it unnatural to have the last 2 be responses. +// It's probably better to modify these functions to take an additional response argument, +// e.g. +// +// template +// bp::object ProjectionEngine::to_map( +// bp::object map, bp::object pbore, bp::object pofs, bp::object det_response, +// bp::object signal, bp::object det_weights, bp::object thread_intervals) +// +// where these would then be extracted using +// +// auto _det_response = BufferWrapper("det_response", det_response, false, vector{n_det,2}); +// +// Could make it optional, but easier to handle this on the python side. +// +// List of functions needing modification: +// * ProjectionEngine::pointing_matrix +// * ProjectionEngine::from_map +// * ProjectionEngine::to_map +// * ProjectionEngine::to_weight_map +// * to_map_single_thread +// * to_weight_map_single_thread +// +// * python/wcs.py/Projectionist +// * python/mapthreads.py/get_threads_domdir + // State the template options class ProjQuat; class ProjFlat; From b790f7c924e281ba15f9fded5654965292680c47 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Sat, 22 Jun 2024 17:10:24 -0400 Subject: [PATCH 02/13] Detector response now supported in c++, and represented in a redone FocalPlane class in python --- demos/demo_proj1.py | 21 +++---- demos/demo_proj2.py | 14 +++-- include/Projection.h | 8 +-- python/proj/coords.py | 126 +++++++++++++++--------------------------- python/proj/wcs.py | 14 ++--- src/Projection.cxx | 72 ++++++++++++++++-------- 6 files changed, 123 insertions(+), 132 deletions(-) diff --git a/demos/demo_proj1.py b/demos/demo_proj1.py index 81337ceb..00b6c188 100644 --- a/demos/demo_proj1.py +++ b/demos/demo_proj1.py @@ -47,6 +47,7 @@ pe = test_utils.get_proj(system, 'TQU', pxz, tiled=args.tiled) ptg = test_utils.get_boresight_quat(system, x, y) ofs = test_utils.get_offsets_quat(system, dx, dy, polphi) +resp= np.ones((n_det,2), np.float32) sig = np.ones((1,n_det,n_t), 'float32') * .5 @@ -119,7 +120,7 @@ def map_delta(a, b): print('Get spin projection factors, too.', end='\n ... ') with Timer() as T: - pix2, spin_proj = pe.pointing_matrix(ptg, ofs, None, None) + pix2, spin_proj = pe.pointing_matrix(ptg, ofs, resp, None, None) del pix, pix_list, pix3, pix2, spin_proj @@ -139,30 +140,30 @@ def map_delta(a, b): else: map1 += np.array([1,0,0])[:,None,None] with Timer() as T: - sig1 = pe.from_map(map1, ptg, ofs, None) + sig1 = pe.from_map(map1, ptg, ofs, resp, None) if 1: print('Project TOD-to-map (TQU)', end='\n ... ') map0 = pe.zeros(3) sig_list = [x for x in sig[0]] with Timer() as T: - map1 = pe.to_map(map0,ptg,ofs,sig_list,None,None) + map1 = pe.to_map(map0,ptg,ofs,resp,sig_list,None,None) if 1: print('TOD-to-map again but with None for input map', end='\n ... ') with Timer() as T: - map1 = pe.to_map(None,ptg,ofs,sig_list,None,None) + map1 = pe.to_map(None,ptg,ofs,resp,sig_list,None,None) if 1: print('Project TOD-to-weights (TQU)', end='\n ... ') map0 = pe.zeros((3, 3)) with Timer() as T: - map2 = pe.to_weight_map(map0,ptg,ofs,None,None) + map2 = pe.to_weight_map(map0,ptg,ofs,resp,None,None) if 1: print('TOD-to-weights again but with None for input map', end='\n ... ') with Timer() as T: - map2 = pe.to_weight_map(None,ptg,ofs,None,None) + map2 = pe.to_weight_map(None,ptg,ofs,resp,None,None) print('Compute thread assignments (OMP prep)... ', end='\n ... ') with Timer(): @@ -171,12 +172,12 @@ def map_delta(a, b): if 1: print('TOD-to-map with OMP (%s): ' % n_omp, end='\n ... ') with Timer() as T: - map1o = pe.to_map(None,ptg,ofs,sig_list,None,threads) + map1o = pe.to_map(None,ptg,ofs,resp,sig_list,None,threads) if 1: print('TOD-to-weights with OMP (%s): ' % n_omp, end='\n ... ') with Timer() as T: - map2o = pe.to_weight_map(None,ptg,ofs,None,threads) + map2o = pe.to_weight_map(None,ptg,ofs,resp,None,threads) print('Checking that OMP and non-OMP forward calcs agree: ', end='\n ... ') assert map_delta(map1, map1o) == 0 @@ -187,7 +188,7 @@ def map_delta(a, b): if 1: print('Cache pointing matrix.', end='\n ...') with Timer() as T: - pix_idx, spin_proj = pe.pointing_matrix(ptg, ofs, None, None) + pix_idx, spin_proj = pe.pointing_matrix(ptg, ofs, resp, None, None) pp = test_utils.get_proj_precomp(args.tiled) print('Map-to-TOD using precomputed pointing matrix', @@ -212,7 +213,7 @@ def map_delta(a, b): print('Checking that it agrees with on-the-fly', end='\n ...') - sig1f = pe.from_map(map1, ptg, ofs, None) + sig1f = pe.from_map(map1, ptg, ofs, resp, None) thresh = map1[0].std() * 1e-6 assert max([np.abs(a - b).max() for a, b in zip(sig1f, sig1p)]) < thresh print('yes') diff --git a/demos/demo_proj2.py b/demos/demo_proj2.py index 839b9dfd..ae23207c 100644 --- a/demos/demo_proj2.py +++ b/demos/demo_proj2.py @@ -82,7 +82,7 @@ def get_pixelizor(emap): ptg = test_utils.get_boresight_quat(system, x*np.pi/180, y*np.pi/180) ofs = test_utils.get_offsets_quat(system, dx, dy, polphi) - +resp= np.ones((n_det,2),np.float32) # # Projection tests. @@ -91,7 +91,7 @@ def get_pixelizor(emap): pe = test_utils.get_proj(system, 'TQU')(pxz) # Project the map into time-domain. -sig0 = pe.from_map(beam, ptg, ofs, None) +sig0 = pe.from_map(beam, ptg, ofs, resp, None) sig0 = np.array(sig0) # Add some noise... @@ -118,10 +118,10 @@ def get_pixelizor(emap): # Then back to map. print('Create timestream...', end='\n ... ') with Timer(): - map1 = pe.to_map(None, ptg, ofs, sig1, det_weights, None) + map1 = pe.to_map(None, ptg, ofs, resp, sig1, det_weights, None) # Get the weight map (matrix). -wmap1 = pe.to_weight_map(None, ptg, ofs, det_weights, None) +wmap1 = pe.to_weight_map(None, ptg, ofs, resp, det_weights, None) wmap1[1,0] = wmap1[0,1] # fill in unpopulated entries... wmap1[2,0] = wmap1[0,2] wmap1[2,1] = wmap1[1,2] @@ -151,6 +151,8 @@ def get_pixelizor(emap): sigd = (sig1[0::2,:] - sig1[1::2,:]) / 2 ofsd = (ofs[::2,...]) +respd= resp[::2] + if not args.uniform_weights: det_weights = det_weights[::2] @@ -160,12 +162,12 @@ def get_pixelizor(emap): # Bin the map again... print('to map...', end='\n ... ') with Timer(): - map1d = pe.to_map(None, ptg, ofsd, sigd, det_weights, None) + map1d = pe.to_map(None, ptg, ofsd, respd, sigd, det_weights, None) # Bin the weights again... print('to weight map...', end='\n ... ') with Timer(): - wmap1d = pe.to_weight_map(None, ptg, ofsd, det_weights, None) + wmap1d = pe.to_weight_map(None, ptg, ofsd, respd, det_weights, None) wmap1d[1,0] = wmap1d[0,1] # fill in unpopulated entries again... diff --git a/include/Projection.h b/include/Projection.h index 8890e095..d5c8c215 100644 --- a/include/Projection.h +++ b/include/Projection.h @@ -48,17 +48,17 @@ class ProjectionEngine { bp::object pixels(bp::object pbore, bp::object pofs, bp::object pixel); vector tile_hits(bp::object pbore, bp::object pofs); bp::object tile_ranges(bp::object pbore, bp::object pofs, bp::object tile_lists); - bp::object pointing_matrix(bp::object pbore, bp::object pofs, + bp::object pointing_matrix(bp::object pbore, bp::object pofs, bp::object response, bp::object pixel, bp::object proj); bp::object zeros(bp::object shape); bp::object pixel_ranges(bp::object pbore, bp::object pofs, bp::object map, int n_domain=-1); bp::object from_map(bp::object map, bp::object pbore, bp::object pofs, - bp::object signal); - bp::object to_map(bp::object map, bp::object pbore, bp::object pofs, + bp::object response, bp::object signal); + bp::object to_map(bp::object map, bp::object pbore, bp::object pofs, bp::object response, bp::object signal, bp::object det_weights, bp::object thread_intervals); bp::object to_weight_map(bp::object map, bp::object pbore, bp::object pofs, - bp::object det_weights, bp::object thread_intervals); + bp::object response, bp::object det_weights, bp::object thread_intervals); int comp_count() const; int index_count() const; diff --git a/python/proj/coords.py b/python/proj/coords.py index a4fea0e6..06d816bc 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -208,93 +208,65 @@ def for_horizon(cls, t, az, el, roll=None, site=None, weather=None): ) return self - def coords(self, det_offsets=None, output=None): + def coords(self, fplane=None, output=None): """Get the celestial coordinates of each detector at each time. Arguments: - det_offset: A dict or list or array of detector offset - tuples. If each tuple has 2 elements or 3 elements, the - typles are interpreted as X,Y[,phi] coordinates in the - conventional way. If 4 elements, it's interpreted as a - quaternion. If this argument is None, then the boresight - pointing is returned. + fplane: A FocalPlane object representing the detector + offsets and responses, or None output: An optional structure for holding the results. For that to work, each element of output must support the buffer protocol. Returns: - If the det_offset was passed in as a dict, then a dict with the same - keys is returned. Otherwise a list is returned. For each - detector, the corresponding returned object is an array with - shape (n_samp, 4) where the four components correspond to - longitude, latitude, cos(psi), sin(psi). - + If fplane is not None, then the result will be + [n_samp,{lon,lat,cos2psi,sin2psi}]. Otherwise it will + be [n_det,n_Samp,{lon,lat,cos2psi,sin2psi}] """ # Get a projector, in CAR. p = so3g.ProjEng_CAR_TQU_NonTiled((1, 1, 1., 1., 1., 1.)) # Pre-process the offsets - collapse = (det_offsets is None) + collapse = (fplane is None) if collapse: - det_offsets = np.array([[1., 0., 0., 0.]]) + fplane = FocalPlane() if output is not None: output = output[None] - redict = isinstance(det_offsets, dict) - if redict: - keys, det_offsets = zip(*det_offsets.items()) - if isinstance(det_offsets[0], quat.quat): - # Individual quat doesn't array() properly... - det_offsets = np.array(quat.G3VectorQuat(det_offsets)) - else: - det_offsets = np.array(det_offsets) - if isinstance(output, dict): - output = [output[k] for k in keys] - if np.shape(det_offsets)[1] == 2: - # XY - x, y = det_offsets[:, 0], det_offsets[:, 1] - theta = np.arcsin((x**2 + y**2)**0.5) - v = np.array([-y, x]) / theta - v[np.isnan(v)] = 0. - det_offsets = np.transpose([np.cos(theta/2), - np.sin(theta/2) * v[0], - np.sin(theta/2) * v[1], - np.zeros(len(theta))]) - det_offsets = quat.G3VectorQuat(det_offsets.copy()) - output = p.coords(self.Q, det_offsets, output) - if redict: - output = OrderedDict(zip(keys, output)) + output = p.coords(self.Q, fplane.quats, fplane.resps, output) if collapse: output = output[0] return output +class FocalPlane: + """This class represents the detector positions and intensity and + polarization responses in the focal plane. -class FocalPlane(OrderedDict): - """This class collects the focal plane positions and orientations for - multiple named detectors. The classmethods can be used to - construct the object from some common input formats. - + This used to be a subclass of OrderedDict, but this was hard to + generalize to per-detector polarization efficiency. """ - @classmethod - def from_xieta(cls, names, xi, eta, gamma=0): - """Creates a FocalPlane object for a set of detector positions - provided in xieta projection plane coordinates. - - Args: - names: vector of detector names. - xi: vector of xi positions, in radians. - eta: vector of eta positions, in radians. - gamma: vector or scalar detector orientation, in radians. - - The (xi,eta) coordinates are Cartesian projection plane - coordinates. When looking at the sky along the telescope - boresight, xi parallel to increasing azimuth and eta is - parallel to increasing elevation. The angle gamma, which - specifies the angle of polarization sensitivity, is measured - from the eta axis, increasing towards the xi axis. - - """ - qs = quat.rotation_xieta(np.asarray(xi), np.asarray(eta), np.asarray(gamma)) - return cls([(n,q) for n, q in zip(names, qs)]) + def __init__(self, quats=None, resps=None): + if quats is None: quats = np.array([[1.0,0.0,0.0,0.0]]) + # Building them this way ensures that + # quats will be an quat coeff array-2 and resps will be a numpy + # array with the right shape, so we don't need to check + # for this when we use FocalPlane later + self.quats = np.array(quat.G3VectorQuat(quats)) + self.resps = np.ones((len(self.quats),2),np.float32) + if resps is not None: + self.resps[:] = resps + @classmethod + def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0): + # The underlying code wants polangle gamma and the T and P + # response, but we support speifying these as the T, Q and U + # response too. Writing it like this handles both cases, and + # as a bonus(?) any combination of them + gamma = gamma + np.arctan2(U,Q)/2 + P = P * (Q**2+U**2)**0.5 + quats = quat.rotation_xieta(np.asarray(xi), np.asarray(eta), np.asarray(gamma)) + resps = np.ones((len(quats),2)) + resps[:,0] = T + resps[:,1] = P + return cls(quats, resps=resps) class Assembly: """This class groups together boresight and detector offset @@ -305,12 +277,11 @@ class Assembly: facilitate more complex arrangements, eventually. """ - def __init__(self, keyed=False, collapse=False): - self.keyed = keyed + def __init__(self, collapse=False): self.collapse = collapse @classmethod - def attach(cls, sight_line, det_offsets): + def attach(cls, sight_line, fplane): """Create an Assembly based on a CelestialSightLine and a FocalPlane. Args: @@ -318,25 +289,16 @@ def attach(cls, sight_line, det_offsets): orientation of the boresight. This can just be a G3VectorQuat if you don't have a whole CelestialSightLine handy. - det_offsets (FocalPlane): The "offsets" of each detector - relative to the boresight. - + fplane (FocalPlane): The "offsets" of each detector + relative to the boresight, and their response to + intensity and polarization """ - keyed = isinstance(det_offsets, dict) self = cls(keyed=keyed) if isinstance(sight_line, quat.G3VectorQuat): self.Q = sight_line else: self.Q = sight_line.Q - if self.keyed: - self.keys = list(det_offsets.keys()) - self.dets = [det_offsets[k] for k in self.keys] - else: - self.dets = det_offsets - # Make sure it's a numpy array. This is dumb. - if isinstance(self.dets[0], quat.quat): - self.dets = quat.G3VectorQuat(self.dets) - self.dets = np.array(self.dets) + self.fplane = fplane return self @classmethod @@ -347,5 +309,5 @@ def for_boresight(cls, sight_line): """ self = cls(collapse=True) self.Q = sight_line.Q - self.dets = [np.array([1., 0, 0, 0])] + self.fplane = FocalPlane() return self diff --git a/python/proj/wcs.py b/python/proj/wcs.py index 49eac62f..f8b3a4df 100644 --- a/python/proj/wcs.py +++ b/python/proj/wcs.py @@ -373,7 +373,7 @@ def get_pointing_matrix(self, assembly): """ projeng = self.get_ProjEng('TQU') q_native = self._cache_q_fp_to_native(assembly.Q) - return projeng.pointing_matrix(q_native, assembly.dets, None, None) + return projeng.pointing_matrix(q_native, assembly.fplane.quats, assembly.fplane.resps, None, None) def get_coords(self, assembly, use_native=False, output=None): """Get the spherical coordinates for the provided pointing Assembly. @@ -446,8 +446,8 @@ def to_map(self, signal, assembly, output=None, det_weights=None, comps = self._guess_comps(output.shape) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - map_out = projeng.to_map( - output, q_native, assembly.dets, signal, det_weights, threads) + map_out = projeng.to_map(output, q_native, assembly.fplane.quats, + assembly.fplane.resps, signal, det_weights, threads) return map_out def to_weights(self, assembly, output=None, det_weights=None, @@ -471,8 +471,8 @@ def to_weights(self, assembly, output=None, det_weights=None, comps = self._guess_comps(output.shape[1:]) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - map_out = projeng.to_weight_map( - output, q_native, assembly.dets, det_weights, threads) + map_out = projeng.to_weight_map(output, q_native, assembly.fplane.quats, + assembly.fplane.resps, det_weights, threads) return map_out def from_map(self, src_map, assembly, signal=None, comps=None): @@ -494,8 +494,8 @@ def from_map(self, src_map, assembly, signal=None, comps=None): comps = self._guess_comps(src_map.shape) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - signal_out = projeng.from_map( - src_map, q_native, assembly.dets, signal) + signal_out = projeng.from_map(src_map, q_native, assembly.fplane.quats, + assembly.fplane.resps, signal) return signal_out def assign_threads(self, assembly, method='domdir', n_threads=None): diff --git a/src/Projection.cxx b/src/Projection.cxx index ae7951de..d5869ae7 100644 --- a/src/Projection.cxx +++ b/src/Projection.cxx @@ -22,6 +22,7 @@ using namespace std; #include "exceptions.h" #include "so_linterp.h" +#include // TRIG_TABLE_SIZE // @@ -855,35 +856,47 @@ inline int Pixelizor2_Flat::GetPixels(int i_det, int i_time, co return iout; } +// This struct is used for representing the total intensity and polarization +// response of a single detector. Originally this was just a length-2 array, +// but those can't easily be returned from functions. If I could ensure +// the buffer had stride-1 in the last dimension, then one could simply pass +// in the address to the first element for a detector, BufferWrapper doesn't make +// that straightforward, so I think this approach is best. +struct Response { FSIGNAL T, P; }; +Response get_response(BufferWrapper & buf, int i_det) { + Response r = {*buf.ptr_2d(i_det,0), *buf.ptr_2d(i_det,1)}; + return r; +} + template static inline -void spin_proj_factors(const double* coords, FSIGNAL *projfacs); +void spin_proj_factors(const double* coords, const Response & response, FSIGNAL *projfacs); template <> -inline void spin_proj_factors(const double* coords, FSIGNAL *projfacs) +inline void spin_proj_factors(const double* coords, const Response & response, FSIGNAL *projfacs) { - projfacs[0] = 1.; + projfacs[0] = response.T; } template <> inline -void spin_proj_factors(const double* coords, FSIGNAL *projfacs) +void spin_proj_factors(const double* coords, const Response & response, FSIGNAL *projfacs) { const double c = coords[2]; const double s = coords[3]; - projfacs[0] = c*c - s*s; - projfacs[1] = 2*c*s; + projfacs[0] = response.P*(c*c - s*s); + projfacs[1] = response.P*(2*c*s); } template <> inline -void spin_proj_factors(const double* coords, FSIGNAL *projfacs) +void spin_proj_factors(const double* coords, const Response & response, FSIGNAL *projfacs) { const double c = coords[2]; const double s = coords[3]; - projfacs[0] = 1.; - projfacs[1] = c*c - s*s; - projfacs[2] = 2*c*s; + projfacs[0] = response.T; + projfacs[1] = response.P*(c*c - s*s); + projfacs[2] = response.P*(2*c*s); } @@ -1068,7 +1081,7 @@ bp::object ProjectionEngine::pixels( // sample template bp::object ProjectionEngine::pointing_matrix( - bp::object pbore, bp::object pofs, bp::object pixel, bp::object proj) + bp::object pbore, bp::object pofs, bp::object response, bp::object pixel, bp::object proj) { auto _none = bp::object(); @@ -1081,9 +1094,10 @@ bp::object ProjectionEngine::pointing_matrix( auto pixel_buf_man = SignalSpace( pixel, "pixel", NPY_INT32, n_det, n_time, P::index_count); - auto proj_buf_man = SignalSpace( proj, "proj", FSIGNAL_NPY_TYPE, n_det, n_time, S::comp_count); + auto _response = BufferWrapper( + "response", response, false, vector{n_det,2}); #pragma omp parallel for for (int i_det = 0; i_det < n_det; ++i_det) { @@ -1093,12 +1107,13 @@ bp::object ProjectionEngine::pointing_matrix( FSIGNAL* const proj_buf = proj_buf_man.data_ptr[i_det]; const int step = pixel_buf_man.steps[0]; int pixel_offset[P::index_count] = {-1}; + Response resp = get_response(_response, i_det); for (int i_time = 0; i_time < n_time; ++i_time) { double coords[4]; FSIGNAL pf[S::comp_count]; pointer.GetCoords(i_det, i_time, (double*)dofs, (double*)coords); _pixelizor.GetPixel(i_det, i_time, (double*)coords, pixel_offset); - spin_proj_factors(coords, pf); + spin_proj_factors(coords, resp, pf); for (int i_dim = 0; i_dim < P::index_count; i_dim++) pix_buf[i_time * pixel_buf_man.steps[0] + @@ -1403,7 +1418,7 @@ bp::object ProjectionEngine::zeros(bp::object shape) template bp::object ProjectionEngine::from_map( - bp::object map, bp::object pbore, bp::object pofs, bp::object signal) + bp::object map, bp::object pbore, bp::object pofs, bp::object response, bp::object signal) { auto _none = bp::object(); @@ -1419,6 +1434,8 @@ bp::object ProjectionEngine::from_map( // Get pointers to the signal and (optional) per-det weights. auto _signalspace = SignalSpace( signal, "signal", FSIGNAL_NPY_TYPE, n_det, n_time); + auto _response = BufferWrapper( + "response", response, false, vector{n_det,2}); #pragma omp parallel for for (int i_det = 0; i_det < n_det; ++i_det) { @@ -1427,10 +1444,11 @@ bp::object ProjectionEngine::from_map( int pixinds[P::interp_count][P::index_count] = {-1}; FSIGNAL pixweights[P::interp_count]; FSIGNAL pf[S::comp_count]; + Response resp = get_response(_response, i_det); for (int i_time = 0; i_time < n_time; ++i_time) { double coords[4]; pointer.GetCoords(i_det, i_time, (double*)dofs, (double*)coords); - spin_proj_factors(coords, pf); + spin_proj_factors(coords, resp, pf); FSIGNAL *sig = (_signalspace.data_ptr[i_det] + _signalspace.steps[0]*i_time); int n_point = _pixelizor.GetPixels(i_det, i_time, coords, pixinds, pixweights); for(int i_point = 0; i_point < n_point; ++i_point) @@ -1445,6 +1463,7 @@ bp::object ProjectionEngine::from_map( template static void to_map_single_thread(Pointer &pointer, + BufferWrapper & _response, P &_pixelizor, const vector & ivals, BufferWrapper &_det_weights, @@ -1462,13 +1481,14 @@ void to_map_single_thread(Pointer &pointer, double coords[4]; FSIGNAL pf[S::comp_count]; pointer.InitPerDet(i_det, dofs); + Response resp = get_response(_response, i_det); // Pointing matrix interpolation stuff int pixinds[P::interp_count][P::index_count] = {-1}; FSIGNAL pixweights[P::interp_count] = {0}; for (auto const &rng: ivals[i_det].segments) { for (int i_time = rng.first; i_time < rng.second; ++i_time) { pointer.GetCoords(i_det, i_time, (double*)dofs, (double*)coords); - spin_proj_factors(coords, pf); + spin_proj_factors(coords, resp, pf); const FSIGNAL sig = _signalspace->data_ptr[i_det][_signalspace->steps[0]*i_time]; // In interpolated mapmaking like bilinear mampamking, each sample hits multipe // pixels, each with its own weight. @@ -1484,6 +1504,7 @@ void to_map_single_thread(Pointer &pointer, template static void to_weight_map_single_thread(Pointer &pointer, + BufferWrapper & _response, P &_pixelizor, vector ivals, BufferWrapper &_det_weights) @@ -1500,12 +1521,13 @@ void to_weight_map_single_thread(Pointer &pointer, double coords[4]; FSIGNAL pf[S::comp_count]; pointer.InitPerDet(i_det, dofs); + Response resp = get_response(_response, i_det); int pixinds[P::interp_count][P::index_count] = {-1}; FSIGNAL pixweights[P::interp_count] = {0}; for (auto const &rng: ivals[i_det].segments) { for (int i_time = rng.first; i_time < rng.second; ++i_time) { pointer.GetCoords(i_det, i_time, (double*)dofs, (double*)coords); - spin_proj_factors(coords, pf); + spin_proj_factors(coords, resp, pf); int n_point = _pixelizor.GetPixels(i_det, i_time, coords, pixinds, pixweights); // Is this enough? Do we need a double loop over i_point? @@ -1592,8 +1614,8 @@ vector>> derive_ranges( template bp::object ProjectionEngine::to_map( - bp::object map, bp::object pbore, bp::object pofs, bp::object signal, bp::object det_weights, - bp::object thread_intervals) + bp::object map, bp::object pbore, bp::object pofs, bp::object response, + bp::object signal, bp::object det_weights, bp::object thread_intervals) { //Initialize it / check inputs. auto pointer = Pointer(); @@ -1613,6 +1635,8 @@ bp::object ProjectionEngine::to_map( signal, "signal", FSIGNAL_NPY_TYPE, n_det, n_time); auto _det_weights = BufferWrapper( "det_weights", det_weights, true, vector{n_det}); + auto _response = BufferWrapper( + "response", response, false, vector{n_det,2}); // For multi-threading, the principle here is that we loop serially // over bunches, and then inside each block all threads loop over @@ -1627,15 +1651,15 @@ bp::object ProjectionEngine::to_map( // or if ivals.size() == 1 #pragma omp parallel for for (int i_thread = 0; i_thread < ivals.size(); i_thread++) - to_map_single_thread(pointer, _pixelizor, ivals[i_thread], _det_weights, &_signalspace); + to_map_single_thread(pointer, _response, _pixelizor, ivals[i_thread], _det_weights, &_signalspace); } return map; } template bp::object ProjectionEngine::to_weight_map( - bp::object map, bp::object pbore, bp::object pofs, bp::object det_weights, - bp::object thread_intervals) + bp::object map, bp::object pbore, bp::object pofs, bp::object response, + bp::object det_weights, bp::object thread_intervals) { auto _none = bp::object(); @@ -1655,6 +1679,8 @@ bp::object ProjectionEngine::to_weight_map( // Get pointer to (optional) per-det weights. auto _det_weights = BufferWrapper( "det_weights", det_weights, true, vector{n_det}); + auto _response = BufferWrapper( + "response", response, false, vector{n_det,2}); // For multi-threading, the principle here is that we loop serially // over bunches, and then inside each block all threads loop over @@ -1669,7 +1695,7 @@ bp::object ProjectionEngine::to_weight_map( // or if ivals.size() == 1 #pragma omp parallel for for (int i_thread = 0; i_thread < ivals.size(); i_thread++) - to_weight_map_single_thread(pointer, _pixelizor, ivals[i_thread], _det_weights); + to_weight_map_single_thread(pointer, _response, _pixelizor, ivals[i_thread], _det_weights); } return map; } From ac3e64a37b22063711a562d9225aa2cd477cb38c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Sun, 23 Jun 2024 13:14:22 -0400 Subject: [PATCH 03/13] Arbitrary response support --- python/proj/coords.py | 22 +++++++++++++++++----- python/proj/mapthreads.py | 24 ++++++++++++------------ python/proj/wcs.py | 20 ++++++++++---------- 3 files changed, 39 insertions(+), 27 deletions(-) diff --git a/python/proj/coords.py b/python/proj/coords.py index 06d816bc..093504cc 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -231,7 +231,7 @@ def coords(self, fplane=None, output=None): fplane = FocalPlane() if output is not None: output = output[None] - output = p.coords(self.Q, fplane.quats, fplane.resps, output) + output = p.coords(self.Q, fplane.quats, output) if collapse: output = output[0] return output @@ -253,20 +253,32 @@ def __init__(self, quats=None, resps=None): self.resps = np.ones((len(self.quats),2),np.float32) if resps is not None: self.resps[:] = resps - + if np.any(~np.isfinite(self.quats)): + raise ValueError("nan/inf values in detector quaternions") + if np.any(~np.isfinite(self.resps)): + raise ValueError("nan/inf values in detector responses") @classmethod - def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0): + def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): # The underlying code wants polangle gamma and the T and P # response, but we support speifying these as the T, Q and U # response too. Writing it like this handles both cases, and # as a bonus(?) any combination of them gamma = gamma + np.arctan2(U,Q)/2 P = P * (Q**2+U**2)**0.5 - quats = quat.rotation_xieta(np.asarray(xi), np.asarray(eta), np.asarray(gamma)) + if hwp: gamma = -gamma + # Broadcast everything to 1d + xi, eta, gamma, T, P, _ = np.broadcast_arrays(xi, eta, gamma, T, P, [0]) + quats = quat.rotation_xieta(xi, eta, gamma) resps = np.ones((len(quats),2)) resps[:,0] = T resps[:,1] = P return cls(quats, resps=resps) + def __repr__(self): + return "FocalPlane(quats=%s,resps=%s)" % (repr(self.quats), repr(self.resps)) + @property + def ndet(self): return len(self.quats) + def __getitem__(self, sel): + return FocalPlane(quats=self.quats[sel], resps=self.resps[sel]) class Assembly: """This class groups together boresight and detector offset @@ -293,7 +305,7 @@ def attach(cls, sight_line, fplane): relative to the boresight, and their response to intensity and polarization """ - self = cls(keyed=keyed) + self = cls() if isinstance(sight_line, quat.G3VectorQuat): self.Q = sight_line else: diff --git a/python/proj/mapthreads.py b/python/proj/mapthreads.py index 0440ea61..92e71f2e 100644 --- a/python/proj/mapthreads.py +++ b/python/proj/mapthreads.py @@ -19,8 +19,8 @@ def get_num_threads(n_threads=None): return so3g.useful_info()['omp_num_threads'] return n_threads -def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None, - active_tiles=None, n_threads=None, offs_rep=None, +def get_threads_domdir(sight, fplane, shape, wcs, tile_shape=None, + active_tiles=None, n_threads=None, fplane_rep=None, plot_prefix=None): """Assign samples to threads according to the dominant scan direction. @@ -37,7 +37,7 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None, Arguments: sight (CelestialSightLine): The boresight pointing. - offs (array of quaternions): The detector pointing. + fplane (FocalPlane): The detector pointing shape (tuple): The map shape. wcs (wcs): The map WCS. tile_shape (tuple): The tile shape, if this should be done using @@ -46,9 +46,9 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None, will be computed along the way. n_threads (int): The number of threads to target (defaults to OMP_NUM_THREADS). - offs_rep (array of quaternions): A representative set of + fplane_rep (FocalPlane): A representative set of detectors, for determining scan direction (if not present then - offs is used for this). + fplane is used for this). plot_prefix (str): If present, pylab will be imported and some plots saved with this name as prefix. @@ -67,8 +67,8 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None, if _m is not None] n_threads = get_num_threads(n_threads) - if offs_rep is None: - offs_rep = offs + if fplane_rep is None: + fplane_rep = fplane if tile_shape is None: # Let's pretend it is, though; this simplifies logic below and @@ -77,7 +77,7 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None, active_tiles = [0] # The full assembly, for later. - asm_full = so3g.proj.Assembly.attach(sight, offs) + asm_full = so3g.proj.Assembly.attach(sight, fplane) # Get a Projectionist -- note it can be used with full or # representative assembly. @@ -93,10 +93,10 @@ def get_threads_domdir(sight, offs, shape, wcs, tile_shape=None, # For the scan direction map, use the "representative" subset # detectors, with polarization direction aligned parallel to # elevation. - xi, eta, gamma = so3g.proj.quat.decompose_xieta(offs_rep) - offs_xl = np.array(so3g.proj.quat.rotation_xieta(xi, eta, gamma*0 + 90*so3g.proj.DEG)) - asm_rep = so3g.proj.Assembly.attach(sight, offs_xl) - sig = np.ones((len(offs_xl), len(asm_rep.Q)), dtype='float32') + xi, eta, gamma = so3g.proj.quat.decompose_xieta(fplane_rep.quats) + fplane_xl = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma*0+90*so3g.proj.DEG) + asm_rep = so3g.proj.Assembly.attach(sight, fplane_xl) + sig = np.ones((fplane_xl.ndet, len(asm_rep.Q)), dtype='float32') scan_maps = pmat.to_map(sig, asm_rep, comps='TQU') # Compute the scan angle, based on Q and U weights. This assumes diff --git a/python/proj/wcs.py b/python/proj/wcs.py index f8b3a4df..4ed06bb9 100644 --- a/python/proj/wcs.py +++ b/python/proj/wcs.py @@ -357,7 +357,7 @@ def get_pixels(self, assembly): """ projeng = self.get_ProjEng('TQU') q_native = self._cache_q_fp_to_native(assembly.Q) - return projeng.pixels(q_native, assembly.dets, None) + return projeng.pixels(q_native, assembly.fplane.quats, None) def get_pointing_matrix(self, assembly): """Get the pointing matrix information, which is to say both the pixel @@ -399,7 +399,7 @@ def get_coords(self, assembly, use_native=False, output=None): q_native = self._cache_q_fp_to_native(assembly.Q) else: q_native = assembly.Q - return projeng.coords(q_native, assembly.dets, output) + return projeng.coords(q_native, assembly.fplane.quats, output) def get_planar(self, assembly, output=None): """Get projection plane coordinates for all detectors at all times. @@ -414,7 +414,7 @@ def get_planar(self, assembly, output=None): """ q_native = self._cache_q_fp_to_native(assembly.Q) projeng = self.get_ProjEng('TQU') - return projeng.coords(q_native, assembly.dets, output) + return projeng.coords(q_native, assembly.fplane.quats, output) def zeros(self, super_shape): """Return a map, filled with zeros, with shape (super_shape,) + @@ -538,11 +538,11 @@ def assign_threads(self, assembly, method='domdir', n_threads=None): if method == 'simple': projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) - omp_ivals = projeng.pixel_ranges(q_native, assembly.dets, None, n_threads) + omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, None, n_threads) return wrap_ivals(omp_ivals) elif method == 'domdir': - offs_rep = assembly.dets[::100] + fplane_rep = assembly.fplane[::100] if (self.tiling is not None) and (self.active_tiles is None): tile_info = self.get_active_tiles(assembly) active_tiles = tile_info['active_tiles'] @@ -550,9 +550,9 @@ def assign_threads(self, assembly, method='domdir', n_threads=None): else: active_tiles = None return mapthreads.get_threads_domdir( - assembly.Q, assembly.dets, shape=self.naxis[::-1], wcs=self.wcs, + assembly.Q, assembly.fplane, shape=self.naxis[::-1], wcs=self.wcs, tile_shape=self.tile_shape, active_tiles=active_tiles, - n_threads=n_threads, offs_rep=offs_rep) + n_threads=n_threads, fplane_rep=fplane_rep) elif method == 'tiles': tile_info = self.get_active_tiles(assembly, assign=n_threads) @@ -578,7 +578,7 @@ def assign_threads_from_map(self, assembly, tmap, n_threads=None): projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) n_threads = mapthreads.get_num_threads(n_threads) - omp_ivals = projeng.pixel_ranges(q_native, assembly.dets, tmap, n_threads) + omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, tmap, n_threads) return wrap_ivals(omp_ivals) def get_active_tiles(self, assembly, assign=False): @@ -617,7 +617,7 @@ def get_active_tiles(self, assembly, assign=False): projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) # This returns a G3VectorInt of length n_tiles giving count of hits per tile. - hits = np.array(projeng.tile_hits(q_native, assembly.dets)) + hits = np.array(projeng.tile_hits(q_native, assembly.fplane.quats)) tiles = np.nonzero(hits)[0] hits = hits[tiles] if assign is True: @@ -631,7 +631,7 @@ def get_active_tiles(self, assembly, assign=False): group_n[idx] += _n group_tiles[idx].append(_t) # Now paint them into Ranges. - R = projeng.tile_ranges(q_native, assembly.dets, group_tiles) + R = projeng.tile_ranges(q_native, assembly.fplane.quats, group_tiles) R = wrap_ivals(R) return { 'active_tiles': list(tiles), From a2eb4dd51d6bb782fe6abffd88b0121f73ac8ef5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Sun, 1 Dec 2024 14:28:20 -0500 Subject: [PATCH 04/13] Passes tests after adding .items(), updating xieta arguments and replacing asm.dets with asm.fplane.quats --- python/proj/coords.py | 3 +++ test/test_proj_astro.py | 11 +++++------ test/test_proj_eng.py | 5 ++--- test/test_proj_eng_hp.py | 4 ++-- 4 files changed, 12 insertions(+), 11 deletions(-) diff --git a/python/proj/coords.py b/python/proj/coords.py index ddb2a6e1..1e113956 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -279,6 +279,9 @@ def __repr__(self): def ndet(self): return len(self.quats) def __getitem__(self, sel): return FocalPlane(quats=self.quats[sel], resps=self.resps[sel]) + def items(self): + for q, resp in zip(quat.G3VectorQuat(self.quats), self.resps): + yield q, resp class Assembly: """This class groups together boresight and detector offset diff --git a/test/test_proj_astro.py b/test/test_proj_astro.py index 35ebe948..3d488f30 100644 --- a/test/test_proj_astro.py +++ b/test/test_proj_astro.py @@ -195,12 +195,10 @@ def test_horizon(self): def test_focalplane(self): # Focal plane - names = ['a', 'b', 'c'] xi = np.array([1., 0., 0.]) * DEG eta = np.array([0., 0., 1.]) * DEG gamma = np.array([45, 0, -45]) * DEG - fp = so3g.proj.FocalPlane.from_xieta(names, xi, eta, gamma) - print(fp) + fp = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma) # Make a CelestialSightLine with some known equatorial pointing. r = so3g.proj.quat.rotation_lonlat(35*DEG, 80*DEG, 15*DEG) @@ -213,13 +211,14 @@ def test_focalplane(self): coords = csl.coords(fp) # Compare to a more direct computation, done here. - for k, roffset in fp.items(): - print('Check detector %s:' % k) + for i, (roffset, resp) in enumerate(fp.items()): + print('Check detector %d:' % i) r1 = r * roffset + print("AAA", type(r1), r1) lon0, lat0, gamma0 = so3g.proj.quat.decompose_lonlat(r1) print(' - inline computation: ', lon0 / DEG, lat0 / DEG, gamma0 / DEG) - lon1, lat1, cosg, sing = coords[k][0] + lon1, lat1, cosg, sing = coords[i][0] gamma1 = np.arctan2(sing, cosg) print(' - proj computation: ', lon1 / DEG, lat1 / DEG, gamma1 / DEG) diff --git a/test/test_proj_eng.py b/test/test_proj_eng.py index b5c190ff..b812d570 100644 --- a/test/test_proj_eng.py +++ b/test/test_proj_eng.py @@ -35,7 +35,7 @@ def get_scan(): def get_basics(clipped=True): t, az, el = get_scan() csl = proj.CelestialSightLine.az_el(t, az, el, weather='vacuum', site='so') - fp = proj.FocalPlane.from_xieta(['a', 'b'], [0., .1*DEG], [0, .1*DEG]) + fp = proj.FocalPlane.from_xieta([0., .1*DEG], [0, .1*DEG]) asm = proj.Assembly.attach(csl, fp) # And a map ... of where? @@ -73,8 +73,7 @@ def test_00_basic(self): det_weights=np.array([0., 0.], dtype='float32'))[0] assert(np.all(m==0)) - # Raise if pointing invalid. - asm.dets[1,2] = np.nan + asm.fplane.quats[1,2] = np.nan with self.assertRaises(ValueError): p.to_map(sig, asm, comps='T') with self.assertRaises(ValueError): diff --git a/test/test_proj_eng_hp.py b/test/test_proj_eng_hp.py index a7eb88c0..6ef20e1a 100644 --- a/test/test_proj_eng_hp.py +++ b/test/test_proj_eng_hp.py @@ -23,7 +23,7 @@ def get_scan(): def get_basics(): t, az, el = get_scan() csl = proj.CelestialSightLine.az_el(t, az, el, weather='vacuum', site='so') - fp = proj.FocalPlane.from_xieta(['a', 'b'], [0., .1*DEG], [0, .1*DEG]) + fp = proj.FocalPlane.from_xieta([0., .1*DEG], [0, .1*DEG]) asm = proj.Assembly.attach(csl, fp) return ((t, az, el), asm) @@ -52,7 +52,7 @@ def test_00_basic(self): assert(np.all(m==0)) # Raise if pointing invalid. - asm.dets[1,2] = np.nan + asm.fplane.quats[1,2] = np.nan with self.assertRaises(ValueError): p.to_map(sig, asm, comps='T') with self.assertRaises(ValueError): From fe20552fc658653cb539ffbd9c72203248725388 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Fri, 10 Jan 2025 08:00:21 -0800 Subject: [PATCH 05/13] Temporary compatibility with old sotodlib --- python/proj/coords.py | 44 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 40 insertions(+), 4 deletions(-) diff --git a/python/proj/coords.py b/python/proj/coords.py index 1e113956..20a23c11 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -243,13 +243,14 @@ class FocalPlane: This used to be a subclass of OrderedDict, but this was hard to generalize to per-detector polarization efficiency. """ - def __init__(self, quats=None, resps=None): - if quats is None: quats = np.array([[1.0,0.0,0.0,0.0]]) + # FIXME: old sotodlib compat - remove dets argument later + def __init__(self, quats=None, resps=None, dets=None): # Building them this way ensures that # quats will be an quat coeff array-2 and resps will be a numpy # array with the right shape, so we don't need to check # for this when we use FocalPlane later - self.quats = np.array(quat.G3VectorQuat(quats)) + if quats is None: self.quats = np.zeros((0,4),np.float64) + else: self.quats = np.array(quat.G3VectorQuat(quats)) self.resps = np.ones((len(self.quats),2),np.float32) if resps is not None: self.resps[:] = resps @@ -257,12 +258,17 @@ def __init__(self, quats=None, resps=None): raise ValueError("nan/inf values in detector quaternions") if np.any(~np.isfinite(self.resps)): raise ValueError("nan/inf values in detector responses") + # FIXME: old sotodlib compat - remove later + self._dets = list(dets) if dets is not None else [] + self._detmap = {name:i for i,name in enumerate(self._dets)} @classmethod def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): # The underlying code wants polangle gamma and the T and P # response, but we support speifying these as the T, Q and U # response too. Writing it like this handles both cases, and # as a bonus(?) any combination of them + # FIXME: old sotodlib compat - remove later + xi, eta, gamma, T, P, Q, U, hwp, dets = cls._xieta_compat(xi,eta,gamma,T,P,Q,U,hwp) gamma = gamma + np.arctan2(U,Q)/2 P = P * (Q**2+U**2)**0.5 if hwp: gamma = -gamma @@ -272,16 +278,43 @@ def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): resps = np.ones((len(quats),2)) resps[:,0] = T resps[:,1] = P - return cls(quats, resps=resps) + # FIXME: old sotodlib compat - remove dets argument later + return cls(quats, resps=resps, dets=dets) def __repr__(self): return "FocalPlane(quats=%s,resps=%s)" % (repr(self.quats), repr(self.resps)) @property def ndet(self): return len(self.quats) + def __len__(self): return self.ndet def __getitem__(self, sel): + # FIXME: old sotodlib compat - remove later + if isinstance(sel, str): return quat.G3VectorQuat(self.quats)[self._dets.index(sel)] return FocalPlane(quats=self.quats[sel], resps=self.resps[sel]) def items(self): for q, resp in zip(quat.G3VectorQuat(self.quats), self.resps): yield q, resp + # FIXME: old sotodlib compat - remove later + @classmethod + def _xieta_compat(cls, xi, eta, gamma, T, P, Q, U, hwp): + # Accept the alternative format (names,xi,eta,gamma) + if isinstance(xi[0], str): + return eta, gamma, T, 1, 1, 1, 0, False, xi + else: + return xi, eta, gamma, T, P, Q, U, hwp, None + # FIXME: old sotodlib compat - remove later + def __setitem__(self, name, q): + # Make coords/pmat.py _get_asm work in old sotodlib. It + # expects to be able to build up a focalplane by assigning + # quats one at a time + q = np.array(quat.G3VectorQuat([q])) + if name in self._detmap: + i = self._detmap[i] + self.quats[i] = q[0] + else: + self._dets.append(name) + self._detmap[name] = len(self._dets)-1 + # This is inefficient, but it's temporary + self.quats = np.concatenate([self.quats,q]) + self.resps = np.concatenate([self.resps,np.ones((1,2),np.float32)]) class Assembly: """This class groups together boresight and detector offset @@ -326,3 +359,6 @@ def for_boresight(cls, sight_line): self.Q = sight_line.Q self.fplane = FocalPlane() return self + # FIXME: old sotodlib compat - remove later + @property + def dets(self): return self.fplane From c396f101be229fb19277a43ec7e96098ab0b5d3d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Fri, 10 Jan 2025 10:29:09 -0800 Subject: [PATCH 06/13] Default FocalPlane is not empty, so make a separate constructor for making one with just the boresight --- python/proj/coords.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/python/proj/coords.py b/python/proj/coords.py index 20a23c11..63d82f9a 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -228,7 +228,7 @@ def coords(self, fplane=None, output=None): # Pre-process the offsets collapse = (fplane is None) if collapse: - fplane = FocalPlane() + fplane = FocalPlane.boresight() if output is not None: output = output[None] output = p.coords(self.Q, fplane.quats, output) @@ -262,6 +262,10 @@ def __init__(self, quats=None, resps=None, dets=None): self._dets = list(dets) if dets is not None else [] self._detmap = {name:i for i,name in enumerate(self._dets)} @classmethod + def boresight(cls): + quats = np.array([[1,0,0,0]],np.float64) + return cls(quats=quats) + @classmethod def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): # The underlying code wants polangle gamma and the T and P # response, but we support speifying these as the T, Q and U @@ -357,7 +361,7 @@ def for_boresight(cls, sight_line): """ self = cls(collapse=True) self.Q = sight_line.Q - self.fplane = FocalPlane() + self.fplane = FocalPlane.boresight() return self # FIXME: old sotodlib compat - remove later @property From 46249922a07473fb63caca30d3e0704c3346e848 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Tue, 14 Jan 2025 05:44:35 -0800 Subject: [PATCH 07/13] Better docstrings and comments --- python/proj/coords.py | 93 ++++++++++++++++++++++++++++++++++++++++--- src/Projection.cxx | 63 +++++++++-------------------- 2 files changed, 105 insertions(+), 51 deletions(-) diff --git a/python/proj/coords.py b/python/proj/coords.py index 63d82f9a..99a392b5 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -214,7 +214,7 @@ def coords(self, fplane=None, output=None): Arguments: fplane: A FocalPlane object representing the detector offsets and responses, or None - output: An optional structure for holding the results. For + output: An optional structure for holding the results. For that to work, each element of output must support the buffer protocol. @@ -240,11 +240,33 @@ class FocalPlane: """This class represents the detector positions and intensity and polarization responses in the focal plane. - This used to be a subclass of OrderedDict, but this was hard to - generalize to per-detector polarization efficiency. + Members: + quats: Array of float64 with shape [ndet,4] representing the + pointing quaternions for each detector + resps: Array of float64 with shape [ndet,2] representing the + total intensity and polarization responsivity of each detector + ndet: The number of detectors (read-only) + + (This used to be a subclass of OrderedDict, but this was hard to + generalize to per-detector polarization efficiency.) """ # FIXME: old sotodlib compat - remove dets argument later def __init__(self, quats=None, resps=None, dets=None): + """Construct a FocalPlane from detector quaternions and responsivities. + + Arguments: + quats: Detector quaternions. Either: + * An array-like of floats with shape [ndet,4] + * An array-like of so3g.proj.quat.quat with shape [ndet] + * An so3g.proj.quat.G3VectorQuat + * None, which results in an empty focalplane with no detectors + resps: Detector responsivities. Either: + * An array-like of floats with shape [ndet,2], where the first and + second number in the last axis are the total intensity and + polarization response respectively + * None, which results in a T and P response of 1 for all detectors. + dets: Deprecated argument temporarily present for backwards + compatibility.""" # Building them this way ensures that # quats will be an quat coeff array-2 and resps will be a numpy # array with the right shape, so we don't need to check @@ -263,10 +285,57 @@ def __init__(self, quats=None, resps=None, dets=None): self._detmap = {name:i for i,name in enumerate(self._dets)} @classmethod def boresight(cls): + """Construct a FocalPlane representing a single detector with a + responsivity of one pointing along the telescope boresight""" quats = np.array([[1,0,0,0]],np.float64) return cls(quats=quats) @classmethod def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): + """Construct a FocalPlane from focal plane coordinates (xi,eta). + + These are Cartesian projection plane coordinates. When looking + at the sky along the telescope boresight, xi is parallel to + increasing azimuth and eta is parallel to increasing elevation. + The angle gamma, which specifies the angle of polarization + sensitivity, is measured from the eta axis, increasing towards + the xi axis. + + Arguments: + xi: An array-like of floats with shape [ndet] + eta: An array-like of floats with shape [ndet] + gamma: The detector polarization angles. [ndet], or a scalar + T, P: The total intensity and polarization responsivity. + Array-like with shape [ndet], or scalar. + Q, U: The Q- and U-polarization responsivity. + Array-like with shape [ndet], or scalar. + Q,U are alternatives to P,gamma for specifying the + polarization responsivity and angle. + + So there are two ways to specify the polarization angle and + responsivity: + 1. gamma and P + 2. Q and U + + Examples, assuming ndet = 2 + * from_xieta(xi, eta, gamma=[0,pi/4]) + Constructs a FocalPlane with T and P responsivity of 1 + and polarization angles of 0 and 45 degrees, representing + a Q-sensitive and U-sensitive detector. + * from_xieta(xi, eta, gamma=[0,pi/4], P=0.5) + Like the above, but with a polarization responsivity of + just 0.5. + * from_xieta(xi, eta, gamma=[0,pi/4], T=[1,0.9], P=[0.5,0.6]) + Like above, but with a detector-dependent intensity and + polarization responsivity. T < P is valid, even though it + wouldn't make sense for most detector types. + * from_xieta(xi, eta, Q=[1,0], U=[0,1]) + Construct the FocalPlane with explicit Q and U responsivity. + This example is equivalent to example 1. + + Usually one would either use gamma,P or Q,U. If they are + combined, then gamma_total = gamma + arctan2(U,Q)/2 and + P_tot = P * (Q**2+U**2)**0.5. + """ # The underlying code wants polangle gamma and the T and P # response, but we support speifying these as the T, Q and U # response too. Writing it like this handles both cases, and @@ -290,10 +359,24 @@ def __repr__(self): def ndet(self): return len(self.quats) def __len__(self): return self.ndet def __getitem__(self, sel): + """Slice the FocalPlane with slice sel, resulting in a new + FocalPlane with a subset of the detectors. Only 1d slice supported, + not integers or multidimensional slices. Example: focal_plane[10:20] + would make a sub-FocalPlane with detector indices 10,11,...,19. + + Deprecated: Temporarily also supports that sel is a detector name, + in which case an spt3g.core.quat is returned for that detector. + This is provided for backwards compatibility.""" # FIXME: old sotodlib compat - remove later if isinstance(sel, str): return quat.G3VectorQuat(self.quats)[self._dets.index(sel)] return FocalPlane(quats=self.quats[sel], resps=self.resps[sel]) def items(self): + """Iterate over detector quaternions and responsivities. Yields + (spt3g.core.quat, array([Tresp,Presp])) pairs. Unlke the raw + entries in the .quats member, which are just numpy arrays with + length 4, spt3g.core.quat are proper quaternon objects that + support quaternion math. + """ for q, resp in zip(quat.G3VectorQuat(self.quats), self.resps): yield q, resp # FIXME: old sotodlib compat - remove later @@ -356,9 +439,7 @@ def attach(cls, sight_line, fplane): @classmethod def for_boresight(cls, sight_line): """Returns an Assembly where a single detector serves as a dummy for - the boresight. - - """ + the boresight.""" self = cls(collapse=True) self.Q = sight_line.Q self.fplane = FocalPlane.boresight() diff --git a/src/Projection.cxx b/src/Projection.cxx index be7e4d7e..96fc7a77 100644 --- a/src/Projection.cxx +++ b/src/Projection.cxx @@ -125,55 +125,28 @@ inline bool isNone(const bp::object &pyo) // Arbitrary detector response support // ----------------------------------- -// Useful to support detectors with polarization response < 1. Also occasionally -// useful with detectors with 0 total intensity response. Currently the detector -// response is only encoded indirectly through the detector angle on the sky, part -// of the per-detector quaternion. Hard to cram the rest of the response into this, -// so must be handled with a new argument. +// It's useful to support detectors with polarization response < 1. It's also +// sometimes useful with detectors with 0 total intensity response. This is +// implemented as follows. // -// spin_proj_factors: cos(2psi),sin(2psi) → detector on-sky TQU response -// psi is here the projected detector angle, which is computed by Pointer::GetCoords -// It makes sense for this to deal with angles since it's the result of a quaternion -// pointing calculation. +// At the top level, the functions where the responsivity is relevant, like +// from_map, to_map etc. take a bp::object reponse argument, representing +// the [ndet,2] detector responsivities. // -// Probably best to leave GetCoords alone, and instead augment spin_proj_factors with -// T and P support. For example: +// This is then converted into a BufferWrapper, before each detector's +// responsivity is extracted using get_response(BufferWrapper & buf, int i_det), +// which returns a Response { FSIGNAL T, P; } struct. // -// void spin_proj_factors(const double * coords, const double * response, FSIGNAL * projfacs) { -// const double c = coords[2]; -// const double s = coords[3]; -// projfacs[0] = response[0]; -// projfacs[1] = response[1]*(c*c - s*s); -// projfacs[2] = response[1]*(2*c*s); -// } +// Finally, this Response struct is passed to spin_proj_factors, e.g. // -// None of the arguments to to_map and its relatives are a natural place to put response. -// It could be shoehorned into pofs by making it length 6, but the first four components -// would hold quaternion coefficients making it unnatural to have the last 2 be responses. -// It's probably better to modify these functions to take an additional response argument, -// e.g. -// -// template -// bp::object ProjectionEngine::to_map( -// bp::object map, bp::object pbore, bp::object pofs, bp::object det_response, -// bp::object signal, bp::object det_weights, bp::object thread_intervals) -// -// where these would then be extracted using -// -// auto _det_response = BufferWrapper("det_response", det_response, false, vector{n_det,2}); -// -// Could make it optional, but easier to handle this on the python side. -// -// List of functions needing modification: -// * ProjectionEngine::pointing_matrix -// * ProjectionEngine::from_map -// * ProjectionEngine::to_map -// * ProjectionEngine::to_weight_map -// * to_map_single_thread -// * to_weight_map_single_thread -// -// * python/wcs.py/Projectionist -// * python/mapthreads.py/get_threads_domdir +// void spin_proj_factors(const double* coords, const Response & response, FSIGNAL *projfacs) +// { +// const double c = coords[2]; +// const double s = coords[3]; +// projfacs[0] = response.T; +// projfacs[1] = response.P*(c*c - s*s); +// projfacs[2] = response.P*(2*c*s); +// } // State the template options class ProjQuat; From bac61f90c19fde75509baa98de67ce5fcf4f2ff8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Tue, 14 Jan 2025 06:53:54 -0800 Subject: [PATCH 08/13] Working on new documentation --- docs/proj.rst | 63 +++++++++++++++++++++++++++++++++++-------- python/proj/coords.py | 5 +++- 2 files changed, 56 insertions(+), 12 deletions(-) diff --git a/docs/proj.rst b/docs/proj.rst index 4ddb42e2..c2b42db8 100644 --- a/docs/proj.rst +++ b/docs/proj.rst @@ -110,26 +110,63 @@ Pointing for many detectors Create a :class:`FocalPlane` object, with some detector positions and orientations:: - names = ['a', 'b', 'c'] - x = np.array([-0.5, 0., 0.5]) * DEG - y = np.zeros(3) + xi = np.array([-0.5, 0.0, 0.5]) * DEG + eta = np.zeros(3) gamma = np.array([0,30,60]) * DEG - fp = so3g.proj.FocalPlane.from_xieta(names, x, y, gamma) + fp = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma) This particular function, :func:`from_xieta() `, will apply the SO standard coordinate -definitions and store a rotation quaternion for each detector. -FocalPlane is just a thinly wrapped OrderedDict, where the detector -name is the key and the value is the rotation quaternion:: +definitions and stores quaternion coefficients (``.quats``) and +responsivities (``.resps``) for each detectors. These are numpy +arrays with shape ``[ndet,4]`` and ``[ndet,2]`` respectively. +``fp.quats[0]`` gives the 4 quaternion coefficients for the +first detector, while ``fp.resps[0]`` gives its total intensity +and polarization responsivity.:: + + >>> fp.quats[2] + array([ 0.86601716, 0.00377878, -0.00218168, 0.49999524]) + >>> fp.resps[2] + array([1., 1.], dtype=float32) + +As you can see, the default responsivity is 1 for both total +intensty and polarization. Note that ``.quats`` contains +quaternion *coefficients*, not quaternion *objects*. To do +quaternion math, you need to convert them to actual quaternion +objects, e.g. ``q = spt3g.core.quat(*fp.quats[0])``, or convert +them all at once with ``qs = spt3g.core.G3VectorQuat(fp.quats)```. + +To represent detectors with responsivity different from 1, +use the ``T`` and ``P`` arguments to :func:`from_xieta()` +to set the total intensity and polarization responsivity +respectively. These can be either single numbers or +array-likes with lengths ``ndet``.:: + + xi = np.array([-0.5, 0.0, 0.5]) * DEG + eta = np.zeros(3) + gamma = np.array([0,30,60]) * DEG + T = np.array([1.0, 0.9, 1.1]) + P = np.array([0.5, 0.6, 0.05]) + fp2 = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma, T=T, P=P) + +Together, gamma, T and P specify the full responsivity of a +detector to the T, Q and U Stokes parameters in focal plane +coordinates. But as an alternative, it's also possible to +specify these directly. The example above is equivalent to:: - >>> fp['c'] - spt3g.core.quat(0.866017,0.00377878,-0.00218168,0.499995) + xi = np.array([-0.5, 0.0, 0.5]) * DEG + eta = np.zeros(3) + gamma = np.array([0,30,60]) * DEG + T = np.array([1.0, 0.9]) + Q = np.array([0.5, 0.3, -0.025]) + U = np.array([0.0, 0.51961524, 0.04330127]) + fp2 = so3g.proj.FocalPlane.from_xieta(xi, eta, T=T, Q=Q, U=U) At this point you could get the celestial coordinates for any one of those detectors:: - # Get vector of quaternion pointings for detector 'a' - q_total = csl.Q * fp['a'] + # Get vector of quaternion pointings for detector 0 + q_total = csl.Q * spt3g.core.quat(*fp.quats[0]) # Decompose q_total into lon, lat, roll angles ra, dec, gamma = so3g.proj.quat.decompose_lonlat(q_total) @@ -144,6 +181,10 @@ to call :func:`coords() ` with the FocalPlane object as first argument:: >>> csl.coords(fp) + [array([[ 0.25715457, -0.92720643, 0.9999161 , 0.01295328]]), + array([[ 0.24261138, -0.92726871, 0.86536489, 0.5011423 ]]), + array([[ 0.22806774, -0.92722945, 0.50890634, 0.86082189]])] + OrderedDict([('a', array([[ 0.22806774, -0.92722945, -0.9999468 , 0.01031487]])), ('b', array([[ 0.24261138, -0.92726871, -0.86536489, -0.5011423 ]])), ('c', array([[ 0.25715457, -0.92720643, -0.48874018, diff --git a/python/proj/coords.py b/python/proj/coords.py index 99a392b5..4cacca1e 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -242,7 +242,10 @@ class FocalPlane: Members: quats: Array of float64 with shape [ndet,4] representing the - pointing quaternions for each detector + pointing quaternions for each detector. Since these are just + plain coefficients, they will need to be converted to quaternion + objects if you want to do quaternion math with them, e.g. + G3VectorQuat(focal_plane.quats). resps: Array of float64 with shape [ndet,2] representing the total intensity and polarization responsivity of each detector ndet: The number of detectors (read-only) From c381023bc1da7fdc28f0f493b90279bea74a4f2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Tue, 14 Jan 2025 07:36:09 -0800 Subject: [PATCH 09/13] More documentation. Updated numbers in examples which seem to have rotted at some point --- docs/proj.rst | 29 ++++++++++++----------------- python/proj/coords.py | 26 ++++++++++++++------------ 2 files changed, 26 insertions(+), 29 deletions(-) diff --git a/docs/proj.rst b/docs/proj.rst index c2b42db8..1d60e373 100644 --- a/docs/proj.rst +++ b/docs/proj.rst @@ -174,7 +174,7 @@ As expected, these coordinates are close to the ones computed before, for the boresight:: >>> print(ra / DEG, dec / DEG) - [13.06731917] [-53.12633433] + [14.73387143] [-53.1250149] But the more expedient way to get pointing for multiple detectors is to call :func:`coords() ` with the @@ -185,11 +185,6 @@ FocalPlane object as first argument:: array([[ 0.24261138, -0.92726871, 0.86536489, 0.5011423 ]]), array([[ 0.22806774, -0.92722945, 0.50890634, 0.86082189]])] - OrderedDict([('a', array([[ 0.22806774, -0.92722945, -0.9999468 , - 0.01031487]])), ('b', array([[ 0.24261138, -0.92726871, -0.86536489, - -0.5011423 ]])), ('c', array([[ 0.25715457, -0.92720643, -0.48874018, - -0.87242939]]))]) - To be clear, ``coords()`` now returns a dictionary whose keys are the detector names. Each value is an array with shape (n_time,4), and at each time step the 4 elements of the array are: ``(lon, lat, @@ -256,17 +251,17 @@ from the pixels in pmap, which are the coordinates of the centers of the pixels:: >>> [x/DEG for x in pix_ra] - [array([13.04], dtype=float32), array([13.88], dtype=float32), array([14.72], dtype=float32)] + [array([14.740001], dtype=float32), array([13.900001], dtype=float32), array([13.059999], dtype=float32)] >>> [x/DEG for x in pix_dec] - [array([-53.100002], dtype=float32), array([-53.100002], dtype=float32), array([-53.100002], dtype=float32)] + [array([-53.12], dtype=float32), array([-53.12], dtype=float32), array([-53.12], dtype=float32)] If you are not getting what you expect, you can grab the pixel indices inferred by the projector -- perhaps your pointing is taking you off the map (in which case the pixel indices would return value -1):: >>> p.get_pixels(asm) - [array([[ 45, 148]], dtype=int32), array([[ 45, 106]], dtype=int32), - array([[45, 64]], dtype=int32)] + [array([[44, 63]], dtype=int32), array([[ 44, 105]], dtype=int32), + array([[ 44, 147]], dtype=int32)] Let's project signal into an intensity map using :func:`Projectionist.to_map`:: @@ -281,9 +276,9 @@ Inspecting the map, we see our signal values occupy the three non-zero pixels: >>> map_out.nonzero() - (array([0, 0, 0]), array([45, 45, 45]), array([ 64, 106, 148])) + (array([0, 0, 0]), array([44, 44, 44]), array([ 63, 105, 147])) >>> map_out[map_out!=0] - array([100., 10., 1.]) + array([ 1., 10., 100.]) If we run this projection again, but pass in this map as a starting @@ -292,7 +287,7 @@ point, the signal will be added to the map: >>> p.to_map(signal, asm, output=map_out, comps='T') array([[[0., 0., 0., ..., 0.]]]) >>> map_out[map_out!=0] - array([200., 20., 2.]) + array([ 2., 20., 200.]) If we instead want to treat the signal as coming from polarization-sensitive detectors, we can request components @@ -305,8 +300,8 @@ according to the projected detector angle on the sky:: >>> map_pol.shape (3, 100, 200) - >>> map_pol[:,45,106] - array([10. , 4.97712803, 8.673419 ]) + >>> map_pol[:,44,105] + array([10., 4.97712803, 8.673419]) For the most basic map-making, the other useful operation is the :func:`Projectionist.to_weights` method. This is used to compute @@ -330,7 +325,7 @@ the upper diagonal has been filled in, for efficiency reasons...:: >>> weight_out.shape (3, 3, 100, 200) - >>> weight_out[...,45,106] + >>> weight_out[...,44,105] array([[1. , 0.49771279, 0.86734194], [0. , 0.24771802, 0.43168718], [0. , 0. , 0.75228202]]) @@ -393,7 +388,7 @@ Inspecting:: >>> threads RangesMatrix(4,3,1) - >>> map_pol2[:,45,106] + >>> map_pol2[:,44,105] array([10. , 4.97712898, 8.67341805]) The same ``threads`` result can be passed to ``p.to_weights``. diff --git a/python/proj/coords.py b/python/proj/coords.py index 4cacca1e..2b5bdd0d 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -320,24 +320,26 @@ def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): 2. Q and U Examples, assuming ndet = 2 - * from_xieta(xi, eta, gamma=[0,pi/4]) + * ``from_xieta(xi, eta, gamma=[0,pi/4])`` Constructs a FocalPlane with T and P responsivity of 1 and polarization angles of 0 and 45 degrees, representing a Q-sensitive and U-sensitive detector. - * from_xieta(xi, eta, gamma=[0,pi/4], P=0.5) + * ``from_xieta(xi, eta, gamma=[0,pi/4], P=0.5)`` Like the above, but with a polarization responsivity of just 0.5. - * from_xieta(xi, eta, gamma=[0,pi/4], T=[1,0.9], P=[0.5,0.6]) + * ``from_xieta(xi, eta, gamma=[0,pi/4], T=[1,0.9], P=[0.5,0.6])`` Like above, but with a detector-dependent intensity and - polarization responsivity. T < P is valid, even though it - wouldn't make sense for most detector types. - * from_xieta(xi, eta, Q=[1,0], U=[0,1]) + polarization responsivity. There is no restriction that + T > P. For the pseudo-detector timestreams one gets after + HWP demodulation, one would have T=0 for the cos-modulated + and sin-modulated timestreams, for example. + * ``from_xieta(xi, eta, Q=[1,0], U=[0,1])`` Construct the FocalPlane with explicit Q and U responsivity. This example is equivalent to example 1. Usually one would either use gamma,P or Q,U. If they are - combined, then gamma_total = gamma + arctan2(U,Q)/2 and - P_tot = P * (Q**2+U**2)**0.5. + combined, then ``gamma_total = gamma + arctan2(U,Q)/2`` and + ``P_tot = P * (Q**2+U**2)**0.5``. """ # The underlying code wants polangle gamma and the T and P # response, but we support speifying these as the T, Q and U @@ -364,20 +366,20 @@ def __len__(self): return self.ndet def __getitem__(self, sel): """Slice the FocalPlane with slice sel, resulting in a new FocalPlane with a subset of the detectors. Only 1d slice supported, - not integers or multidimensional slices. Example: focal_plane[10:20] + not integers or multidimensional slices. Example: ``focal_plane[10:20]`` would make a sub-FocalPlane with detector indices 10,11,...,19. Deprecated: Temporarily also supports that sel is a detector name, - in which case an spt3g.core.quat is returned for that detector. + in which case an ``spt3g.core.quat`` is returned for that detector. This is provided for backwards compatibility.""" # FIXME: old sotodlib compat - remove later if isinstance(sel, str): return quat.G3VectorQuat(self.quats)[self._dets.index(sel)] return FocalPlane(quats=self.quats[sel], resps=self.resps[sel]) def items(self): """Iterate over detector quaternions and responsivities. Yields - (spt3g.core.quat, array([Tresp,Presp])) pairs. Unlke the raw + ``(spt3g.core.quat, array([Tresp,Presp]))`` pairs. Unlke the raw entries in the .quats member, which are just numpy arrays with - length 4, spt3g.core.quat are proper quaternon objects that + length 4, ``spt3g.core.quat`` are proper quaternon objects that support quaternion math. """ for q, resp in zip(quat.G3VectorQuat(self.quats), self.resps): From 63f92f181e578f2c7e36d46f9833b4a3f4dce13e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Thu, 16 Jan 2025 07:19:24 -0800 Subject: [PATCH 10/13] Matthews comments, plus go to G3VectorQuat instead of numpy array --- docs/proj.rst | 25 +++++++------- python/proj/coords.py | 71 +++++++++++++++++++++++----------------- python/proj/wcs.py | 24 +++++++------- test/test_proj_eng.py | 4 ++- test/test_proj_eng_hp.py | 11 ++++--- 5 files changed, 74 insertions(+), 61 deletions(-) diff --git a/docs/proj.rst b/docs/proj.rst index 1d60e373..6175da34 100644 --- a/docs/proj.rst +++ b/docs/proj.rst @@ -85,9 +85,9 @@ quaternions into rotation angles:: >>> csl.Q spt3g.core.G3VectorQuat([(-0.0384748,0.941783,0.114177,0.313891)]) - >>> csl.coords() + >>> csl.coords() array([[ 0.24261138, -0.92726871, -0.99999913, -0.00131952]]) - + The :func:`coords() ` returns an array with shape (n_time, 4); each 4-tuple contains values ``(lon, lat, cos(gamma), sin(gamma))``. The ``lon`` and ``lat`` are the celestial @@ -117,12 +117,12 @@ orientations:: This particular function, :func:`from_xieta() `, will apply the SO standard coordinate -definitions and stores quaternion coefficients (``.quats``) and -responsivities (``.resps``) for each detectors. These are numpy -arrays with shape ``[ndet,4]`` and ``[ndet,2]`` respectively. -``fp.quats[0]`` gives the 4 quaternion coefficients for the -first detector, while ``fp.resps[0]`` gives its total intensity -and polarization responsivity.:: +definitions and stores quaternions (``.quats``) and +responsivities (``.resps``) for each detector. These are a +``G3VectorQuat`` wth length ``ndet`` and a numpy array +with shape ``[ndet,2]`` respectively. ``fp.quats[0]`` gives +the quaternion for the first detector, while ``fp.resps[0]`` +gives its total intensity and polarization responsivity.:: >>> fp.quats[2] array([ 0.86601716, 0.00377878, -0.00218168, 0.49999524]) @@ -133,7 +133,7 @@ As you can see, the default responsivity is 1 for both total intensty and polarization. Note that ``.quats`` contains quaternion *coefficients*, not quaternion *objects*. To do quaternion math, you need to convert them to actual quaternion -objects, e.g. ``q = spt3g.core.quat(*fp.quats[0])``, or convert +objects, e.g. ``q = fp.quats[0]``, or convert them all at once with ``qs = spt3g.core.G3VectorQuat(fp.quats)```. To represent detectors with responsivity different from 1, @@ -166,7 +166,7 @@ At this point you could get the celestial coordinates for any one of those detectors:: # Get vector of quaternion pointings for detector 0 - q_total = csl.Q * spt3g.core.quat(*fp.quats[0]) + q_total = csl.Q * fp.quats[0] # Decompose q_total into lon, lat, roll angles ra, dec, gamma = so3g.proj.quat.decompose_lonlat(q_total) @@ -185,9 +185,8 @@ FocalPlane object as first argument:: array([[ 0.24261138, -0.92726871, 0.86536489, 0.5011423 ]]), array([[ 0.22806774, -0.92722945, 0.50890634, 0.86082189]])] -To be clear, ``coords()`` now returns a dictionary whose keys are the -detector names. Each value is an array with shape (n_time,4), and at -each time step the 4 elements of the array are: ``(lon, lat, +So ``coords()`` returns a list of numpy arrays with shape (n_time,4), +and at each time step the 4 elements of the array are: ``(lon, lat, cos(gamma), sin(gamma))``. diff --git a/python/proj/coords.py b/python/proj/coords.py index 2b5bdd0d..a2eea6d7 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -219,9 +219,9 @@ def coords(self, fplane=None, output=None): buffer protocol. Returns: - If fplane is not None, then the result will be + If fplane is None, then the result will be [n_samp,{lon,lat,cos2psi,sin2psi}]. Otherwise it will - be [n_det,n_Samp,{lon,lat,cos2psi,sin2psi}] + be [n_det,n_samp,{lon,lat,cos2psi,sin2psi}] """ # Get a projector, in CAR. p = so3g.ProjEng_CAR_TQU_NonTiled((1, 1, 1., 1., 1., 1.)) @@ -231,7 +231,7 @@ def coords(self, fplane=None, output=None): fplane = FocalPlane.boresight() if output is not None: output = output[None] - output = p.coords(self.Q, fplane.quats, output) + output = p.coords(self.Q, fplane.coeffs(), output) if collapse: output = output[0] return output @@ -241,11 +241,9 @@ class FocalPlane: polarization responses in the focal plane. Members: - quats: Array of float64 with shape [ndet,4] representing the - pointing quaternions for each detector. Since these are just - plain coefficients, they will need to be converted to quaternion - objects if you want to do quaternion math with them, e.g. - G3VectorQuat(focal_plane.quats). + quats: G3VectorQuat representing the pointing quaternions for + each detector. Can be turned into a numpy array of coefficients + with np.array(). Or call .coeffs() to get them directly. resps: Array of float64 with shape [ndet,2] representing the total intensity and polarization responsivity of each detector ndet: The number of detectors (read-only) @@ -274,8 +272,10 @@ def __init__(self, quats=None, resps=None, dets=None): # quats will be an quat coeff array-2 and resps will be a numpy # array with the right shape, so we don't need to check # for this when we use FocalPlane later - if quats is None: self.quats = np.zeros((0,4),np.float64) - else: self.quats = np.array(quat.G3VectorQuat(quats)) + if quats is None: quats = [] + # Asarray needed because G3VectorQuat doesn't handle list of lists, + # which we want to be able to accept + self.quats = quat.G3VectorQuat(np.asarray(quats)) self.resps = np.ones((len(self.quats),2),np.float32) if resps is not None: self.resps[:] = resps @@ -286,15 +286,22 @@ def __init__(self, quats=None, resps=None, dets=None): # FIXME: old sotodlib compat - remove later self._dets = list(dets) if dets is not None else [] self._detmap = {name:i for i,name in enumerate(self._dets)} + def coeffs(self): + """Return an [ndet,4] numpy array representing the quaternion + coefficients of ``.quats``. Useful for passing the detector + pointing to C++ code""" + return np.array(self.quats, dtype=np.float64) @classmethod def boresight(cls): """Construct a FocalPlane representing a single detector with a responsivity of one pointing along the telescope boresight""" - quats = np.array([[1,0,0,0]],np.float64) - return cls(quats=quats) + return cls(quats=[[1,0,0,0]]) + # FIXME: old sotodlib compat - expand to actual argument list later @classmethod - def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): - """Construct a FocalPlane from focal plane coordinates (xi,eta). + def from_xieta(cls, *args, **kwargs): + """ + def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): + Construct a FocalPlane from focal plane coordinates (xi,eta). These are Cartesian projection plane coordinates. When looking at the sky along the telescope boresight, xi is parallel to @@ -346,7 +353,7 @@ def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): # response too. Writing it like this handles both cases, and # as a bonus(?) any combination of them # FIXME: old sotodlib compat - remove later - xi, eta, gamma, T, P, Q, U, hwp, dets = cls._xieta_compat(xi,eta,gamma,T,P,Q,U,hwp) + xi, eta, gamma, T, P, Q, U, hwp, dets = cls._xieta_compat(*args, **kwargs) gamma = gamma + np.arctan2(U,Q)/2 P = P * (Q**2+U**2)**0.5 if hwp: gamma = -gamma @@ -359,22 +366,25 @@ def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): # FIXME: old sotodlib compat - remove dets argument later return cls(quats, resps=resps, dets=dets) def __repr__(self): - return "FocalPlane(quats=%s,resps=%s)" % (repr(self.quats), repr(self.resps)) + return "FocalPlane(quats=%s,resps=%s)" % (repr(self.coeffs()), repr(self.resps)) @property def ndet(self): return len(self.quats) def __len__(self): return self.ndet def __getitem__(self, sel): """Slice the FocalPlane with slice sel, resulting in a new - FocalPlane with a subset of the detectors. Only 1d slice supported, - not integers or multidimensional slices. Example: ``focal_plane[10:20]`` - would make a sub-FocalPlane with detector indices 10,11,...,19. + FocalPlane with a subset of the detectors. Only a 1d slice + or boolean mask is supported, not integers or multidimensional + slices. Example: ``focal_plane[10:20]`` would make a + sub-FocalPlane with detector indices 10,11,...,19. Deprecated: Temporarily also supports that sel is a detector name, in which case an ``spt3g.core.quat`` is returned for that detector. This is provided for backwards compatibility.""" # FIXME: old sotodlib compat - remove later - if isinstance(sel, str): return quat.G3VectorQuat(self.quats)[self._dets.index(sel)] - return FocalPlane(quats=self.quats[sel], resps=self.resps[sel]) + if isinstance(sel, str): return self.quats[self._dets.index(sel)] + # We go via .coeffs() here to get around G3VectorQuat's lack + # of boolean mask support + return FocalPlane(quats=self.coeffs()[sel], resps=self.resps[sel]) def items(self): """Iterate over detector quaternions and responsivities. Yields ``(spt3g.core.quat, array([Tresp,Presp]))`` pairs. Unlke the raw @@ -382,30 +392,31 @@ def items(self): length 4, ``spt3g.core.quat`` are proper quaternon objects that support quaternion math. """ - for q, resp in zip(quat.G3VectorQuat(self.quats), self.resps): + for q, resp in zip(self.quats, self.resps): yield q, resp # FIXME: old sotodlib compat - remove later - @classmethod - def _xieta_compat(cls, xi, eta, gamma, T, P, Q, U, hwp): + @staticmethod + def _xieta_compat(*args, **kwargs): # Accept the alternative format (names,xi,eta,gamma) - if isinstance(xi[0], str): - return eta, gamma, T, 1, 1, 1, 0, False, xi + def helper(xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False, dets=None): + return xi, eta, gamma, T, P, Q, U, hwp, dets + if isinstance(args[0][0], str): + return helper(*args[1:], dets=args[0], **kwargs) else: - return xi, eta, gamma, T, P, Q, U, hwp, None + return helper(*args, **kwargs) # FIXME: old sotodlib compat - remove later def __setitem__(self, name, q): # Make coords/pmat.py _get_asm work in old sotodlib. It # expects to be able to build up a focalplane by assigning # quats one at a time - q = np.array(quat.G3VectorQuat([q])) if name in self._detmap: i = self._detmap[i] - self.quats[i] = q[0] + self.quats[i] = q else: self._dets.append(name) self._detmap[name] = len(self._dets)-1 + self.quats.append(q) # This is inefficient, but it's temporary - self.quats = np.concatenate([self.quats,q]) self.resps = np.concatenate([self.resps,np.ones((1,2),np.float32)]) class Assembly: diff --git a/python/proj/wcs.py b/python/proj/wcs.py index 860242e4..74f2753f 100644 --- a/python/proj/wcs.py +++ b/python/proj/wcs.py @@ -183,7 +183,7 @@ def get_pixels(self, assembly): """ projeng = self.get_ProjEng('TQU') q_native = self._cache_q_fp_to_native(assembly.Q) - return projeng.pixels(q_native, assembly.fplane.quats, None) + return projeng.pixels(q_native, assembly.fplane.coeffs(), None) def get_pointing_matrix(self, assembly): """Get the pointing matrix information, which is to say both the pixel @@ -199,7 +199,7 @@ def get_pointing_matrix(self, assembly): """ projeng = self.get_ProjEng('TQU') q_native = self._cache_q_fp_to_native(assembly.Q) - return projeng.pointing_matrix(q_native, assembly.fplane.quats, assembly.fplane.resps, None, None) + return projeng.pointing_matrix(q_native, assembly.fplane.coeffs(), assembly.fplane.resps, None, None) def get_coords(self, assembly, use_native=False, output=None): """Get the spherical coordinates for the provided pointing Assembly. @@ -225,7 +225,7 @@ def get_coords(self, assembly, use_native=False, output=None): q_native = self._cache_q_fp_to_native(assembly.Q) else: q_native = assembly.Q - return projeng.coords(q_native, assembly.fplane.quats, output) + return projeng.coords(q_native, assembly.fplane.coeffs(), output) def get_planar(self, assembly, output=None): """Get projection plane coordinates for all detectors at all times. @@ -240,7 +240,7 @@ def get_planar(self, assembly, output=None): """ q_native = self._cache_q_fp_to_native(assembly.Q) projeng = self.get_ProjEng('TQU') - return projeng.coords(q_native, assembly.fplane.quats, output) + return projeng.coords(q_native, assembly.fplane.coeffs(), output) def zeros(self, super_shape): """Return a map, filled with zeros, with shape (super_shape,) + @@ -272,7 +272,7 @@ def to_map(self, signal, assembly, output=None, det_weights=None, comps = self._guess_comps(self._get_map_shape(output)) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - map_out = projeng.to_map(output, q_native, assembly.fplane.quats, + map_out = projeng.to_map(output, q_native, assembly.fplane.coeffs(), assembly.fplane.resps, signal, det_weights, threads) return map_out @@ -298,7 +298,7 @@ def to_weights(self, assembly, output=None, det_weights=None, comps = self._guess_comps(shape[1:]) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - map_out = projeng.to_weight_map(output, q_native, assembly.fplane.quats, + map_out = projeng.to_weight_map(output, q_native, assembly.fplane.coeffs(), assembly.fplane.resps, det_weights, threads) return map_out @@ -318,7 +318,7 @@ def from_map(self, src_map, assembly, signal=None, comps=None): comps = self._guess_comps(self._get_map_shape(src_map)) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - signal_out = projeng.from_map(src_map, q_native, assembly.fplane.quats, + signal_out = projeng.from_map(src_map, q_native, assembly.fplane.coeffs(), assembly.fplane.resps, signal) return signal_out @@ -362,7 +362,7 @@ def assign_threads(self, assembly, method='domdir', n_threads=None): if method == 'simple': projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) - omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, None, n_threads) + omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.coeffs(), None, n_threads) return wrap_ivals(omp_ivals) elif method == 'domdir': @@ -402,7 +402,7 @@ def assign_threads_from_map(self, assembly, tmap, n_threads=None): projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) n_threads = mapthreads.get_num_threads(n_threads) - omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, tmap, n_threads) + omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.coeffs(), tmap, n_threads) return wrap_ivals(omp_ivals) def get_active_tiles(self, assembly, assign=False): @@ -441,7 +441,7 @@ def get_active_tiles(self, assembly, assign=False): projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) # This returns a G3VectorInt of length n_tiles giving count of hits per tile. - hits = np.array(projeng.tile_hits(q_native, assembly.fplane.quats)) + hits = np.array(projeng.tile_hits(q_native, assembly.fplane.coeffs())) tiles = np.nonzero(hits)[0] hits = hits[tiles] if assign is True: @@ -460,7 +460,7 @@ def get_active_tiles(self, assembly, assign=False): # print(f"Warning: Threads poorly balanced. Max/mean hits = {max_ratio}") # Now paint them into Ranges. - R = projeng.tile_ranges(q_native, assembly.fplane.quats, group_tiles) + R = projeng.tile_ranges(q_native, assembly.fplane.coeffs(), group_tiles) R = wrap_ivals(R) return { 'active_tiles': list(tiles), @@ -795,7 +795,7 @@ def get_coords(self, assembly, use_native=False, output=None): q_native = self._cache_q_fp_to_native(assembly.Q) else: q_native = assembly.Q - return projeng.coords(q_native, assembly.dets, output) + return projeng.coords(q_native, assembly.fplane.coeffs(), output) def from_map(self, src_map, assembly, signal=None, comps=None): """De-project from a map, returning a Signal-like object. diff --git a/test/test_proj_eng.py b/test/test_proj_eng.py index b812d570..17bd261c 100644 --- a/test/test_proj_eng.py +++ b/test/test_proj_eng.py @@ -73,7 +73,9 @@ def test_00_basic(self): det_weights=np.array([0., 0.], dtype='float32'))[0] assert(np.all(m==0)) - asm.fplane.quats[1,2] = np.nan + # Can't assign to quat fields, so do + # it this way instead + asm.fplane.quats[1] = asm.fplane.quats[1]*np.nan with self.assertRaises(ValueError): p.to_map(sig, asm, comps='T') with self.assertRaises(ValueError): diff --git a/test/test_proj_eng_hp.py b/test/test_proj_eng_hp.py index 6ef20e1a..b407ade8 100644 --- a/test/test_proj_eng_hp.py +++ b/test/test_proj_eng_hp.py @@ -32,7 +32,7 @@ class TestProjEngHP(unittest.TestCase): """Test ProjectionistHealpix Based on TestProjEng """ - + def test_00_basic(self): scan, asm = get_basics() nside = 128 @@ -51,8 +51,9 @@ def test_00_basic(self): det_weights=np.array([0., 0.], dtype='float32'))[0] assert(np.all(m==0)) - # Raise if pointing invalid. - asm.fplane.quats[1,2] = np.nan + # Can't assign to quat fields, so do + # it this way instead + asm.fplane.quats[1] = asm.fplane.quats[1]*np.nan with self.assertRaises(ValueError): p.to_map(sig, asm, comps='T') with self.assertRaises(ValueError): @@ -74,7 +75,7 @@ def test_10_tiled(self): assert(np.any(w2)) # Identify active subtiles? print(p.active_tiles) - + def test_20_threads(self): for (tiled, interpol, method) in itertools.product( [False, True], @@ -88,7 +89,7 @@ def test_20_threads(self): nside_tile = 8 else: nside_tile = None - + p = proj.ProjectionistHealpix.for_healpix(nside, nside_tile, interpol=interpol) sig = np.ones((2, len(scan[0])), 'float32') n_threads = 3 From 5efc1f71d50a1116623681ee44069ea986f253c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Thu, 16 Jan 2025 09:23:31 -0800 Subject: [PATCH 11/13] Ready for new review --- demos/demo_proj_hp.py | 11 ++--------- docs/proj.rst | 20 +++++++------------- python/proj/coords.py | 2 +- python/proj/wcs.py | 24 ++++++++++++------------ 4 files changed, 22 insertions(+), 35 deletions(-) diff --git a/demos/demo_proj_hp.py b/demos/demo_proj_hp.py index 2b095b0c..9f81fcb4 100644 --- a/demos/demo_proj_hp.py +++ b/demos/demo_proj_hp.py @@ -191,13 +191,6 @@ def simulate_detectors(ndets, rmax, seed=0): eta = rr * np.sin(phi) return xi, eta, gamma -def make_fp(xi, eta, gamma): - fp = so3g.proj.quat.rotation_xieta(xi, eta, gamma) - so3g_fp = so3g.proj.FocalPlane() - for i, q in enumerate(fp): - so3g_fp[f'a{i}'] = q - return so3g_fp - def extract_coord(sight, fp, isignal=0, groupsize=100): coord_out = [] for ii in range(0, len(fp), groupsize): @@ -208,14 +201,14 @@ def extract_coord(sight, fp, isignal=0, groupsize=100): def make_signal(ndets, rmax, sight, isignal, seed=0): xi, eta, gamma = simulate_detectors(ndets, rmax, seed) - fp = so3g.proj.quat.rotation_xieta(xi, eta, gamma) + fp = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma) signal = extract_coord(sight, fp, isignal) return signal def make_tod(time, az, el, ndets, rmax, isignal, seed=0): sight = so3g.proj.CelestialSightLine.naive_az_el(time, az, el) signal = make_signal(ndets, rmax, sight, isignal, seed) - so3g_fp = make_fp(*simulate_detectors(ndets, rmax, seed)) + so3g_fp = so3g.proj.FocalPlane.from_xieta(*simulate_detectors(ndets, rmax, seed)) asm = so3g.proj.Assembly.attach(sight, so3g_fp) return signal, asm diff --git a/docs/proj.rst b/docs/proj.rst index 6175da34..9401ada6 100644 --- a/docs/proj.rst +++ b/docs/proj.rst @@ -125,22 +125,16 @@ the quaternion for the first detector, while ``fp.resps[0]`` gives its total intensity and polarization responsivity.:: >>> fp.quats[2] - array([ 0.86601716, 0.00377878, -0.00218168, 0.49999524]) + spt3g.core.quat(0.866017,0.00377878,-0.00218168,0.499995) >>> fp.resps[2] array([1., 1.], dtype=float32) As you can see, the default responsivity is 1 for both total -intensty and polarization. Note that ``.quats`` contains -quaternion *coefficients*, not quaternion *objects*. To do -quaternion math, you need to convert them to actual quaternion -objects, e.g. ``q = fp.quats[0]``, or convert -them all at once with ``qs = spt3g.core.G3VectorQuat(fp.quats)```. - -To represent detectors with responsivity different from 1, -use the ``T`` and ``P`` arguments to :func:`from_xieta()` -to set the total intensity and polarization responsivity -respectively. These can be either single numbers or -array-likes with lengths ``ndet``.:: +intensty and polarization. To represent detectors with +responsivity different from 1, use the ``T`` and ``P`` +arguments to :func:`from_xieta()` to set the total intensity +and polarization responsivity respectively. These can be +either single numbers or array-likes with lengths ``ndet``.:: xi = np.array([-0.5, 0.0, 0.5]) * DEG eta = np.zeros(3) @@ -157,7 +151,7 @@ specify these directly. The example above is equivalent to:: xi = np.array([-0.5, 0.0, 0.5]) * DEG eta = np.zeros(3) gamma = np.array([0,30,60]) * DEG - T = np.array([1.0, 0.9]) + T = np.array([1.0, 0.9, 1.1]) Q = np.array([0.5, 0.3, -0.025]) U = np.array([0.0, 0.51961524, 0.04330127]) fp2 = so3g.proj.FocalPlane.from_xieta(xi, eta, T=T, Q=Q, U=U) diff --git a/python/proj/coords.py b/python/proj/coords.py index a2eea6d7..1fcb2ae1 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -231,7 +231,7 @@ def coords(self, fplane=None, output=None): fplane = FocalPlane.boresight() if output is not None: output = output[None] - output = p.coords(self.Q, fplane.coeffs(), output) + output = p.coords(self.Q, fplane.quats), output) if collapse: output = output[0] return output diff --git a/python/proj/wcs.py b/python/proj/wcs.py index 74f2753f..8c2005d3 100644 --- a/python/proj/wcs.py +++ b/python/proj/wcs.py @@ -183,7 +183,7 @@ def get_pixels(self, assembly): """ projeng = self.get_ProjEng('TQU') q_native = self._cache_q_fp_to_native(assembly.Q) - return projeng.pixels(q_native, assembly.fplane.coeffs(), None) + return projeng.pixels(q_native, assembly.fplane.quats, None) def get_pointing_matrix(self, assembly): """Get the pointing matrix information, which is to say both the pixel @@ -199,7 +199,7 @@ def get_pointing_matrix(self, assembly): """ projeng = self.get_ProjEng('TQU') q_native = self._cache_q_fp_to_native(assembly.Q) - return projeng.pointing_matrix(q_native, assembly.fplane.coeffs(), assembly.fplane.resps, None, None) + return projeng.pointing_matrix(q_native, assembly.fplane.quats, assembly.fplane.resps, None, None) def get_coords(self, assembly, use_native=False, output=None): """Get the spherical coordinates for the provided pointing Assembly. @@ -225,7 +225,7 @@ def get_coords(self, assembly, use_native=False, output=None): q_native = self._cache_q_fp_to_native(assembly.Q) else: q_native = assembly.Q - return projeng.coords(q_native, assembly.fplane.coeffs(), output) + return projeng.coords(q_native, assembly.fplane.quats, output) def get_planar(self, assembly, output=None): """Get projection plane coordinates for all detectors at all times. @@ -240,7 +240,7 @@ def get_planar(self, assembly, output=None): """ q_native = self._cache_q_fp_to_native(assembly.Q) projeng = self.get_ProjEng('TQU') - return projeng.coords(q_native, assembly.fplane.coeffs(), output) + return projeng.coords(q_native, assembly.fplane.quats, output) def zeros(self, super_shape): """Return a map, filled with zeros, with shape (super_shape,) + @@ -272,7 +272,7 @@ def to_map(self, signal, assembly, output=None, det_weights=None, comps = self._guess_comps(self._get_map_shape(output)) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - map_out = projeng.to_map(output, q_native, assembly.fplane.coeffs(), + map_out = projeng.to_map(output, q_native, assembly.fplane.quats, assembly.fplane.resps, signal, det_weights, threads) return map_out @@ -298,7 +298,7 @@ def to_weights(self, assembly, output=None, det_weights=None, comps = self._guess_comps(shape[1:]) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - map_out = projeng.to_weight_map(output, q_native, assembly.fplane.coeffs(), + map_out = projeng.to_weight_map(output, q_native, assembly.fplane.quats, assembly.fplane.resps, det_weights, threads) return map_out @@ -318,7 +318,7 @@ def from_map(self, src_map, assembly, signal=None, comps=None): comps = self._guess_comps(self._get_map_shape(src_map)) projeng = self.get_ProjEng(comps) q_native = self._cache_q_fp_to_native(assembly.Q) - signal_out = projeng.from_map(src_map, q_native, assembly.fplane.coeffs(), + signal_out = projeng.from_map(src_map, q_native, assembly.fplane.quats, assembly.fplane.resps, signal) return signal_out @@ -362,7 +362,7 @@ def assign_threads(self, assembly, method='domdir', n_threads=None): if method == 'simple': projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) - omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.coeffs(), None, n_threads) + omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, None, n_threads) return wrap_ivals(omp_ivals) elif method == 'domdir': @@ -402,7 +402,7 @@ def assign_threads_from_map(self, assembly, tmap, n_threads=None): projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) n_threads = mapthreads.get_num_threads(n_threads) - omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.coeffs(), tmap, n_threads) + omp_ivals = projeng.pixel_ranges(q_native, assembly.fplane.quats, tmap, n_threads) return wrap_ivals(omp_ivals) def get_active_tiles(self, assembly, assign=False): @@ -441,7 +441,7 @@ def get_active_tiles(self, assembly, assign=False): projeng = self.get_ProjEng('T') q_native = self._cache_q_fp_to_native(assembly.Q) # This returns a G3VectorInt of length n_tiles giving count of hits per tile. - hits = np.array(projeng.tile_hits(q_native, assembly.fplane.coeffs())) + hits = np.array(projeng.tile_hits(q_native, assembly.fplane.quats)) tiles = np.nonzero(hits)[0] hits = hits[tiles] if assign is True: @@ -460,7 +460,7 @@ def get_active_tiles(self, assembly, assign=False): # print(f"Warning: Threads poorly balanced. Max/mean hits = {max_ratio}") # Now paint them into Ranges. - R = projeng.tile_ranges(q_native, assembly.fplane.coeffs(), group_tiles) + R = projeng.tile_ranges(q_native, assembly.fplane.quats, group_tiles) R = wrap_ivals(R) return { 'active_tiles': list(tiles), @@ -795,7 +795,7 @@ def get_coords(self, assembly, use_native=False, output=None): q_native = self._cache_q_fp_to_native(assembly.Q) else: q_native = assembly.Q - return projeng.coords(q_native, assembly.fplane.coeffs(), output) + return projeng.coords(q_native, assembly.fplane.quats, output) def from_map(self, src_map, assembly, signal=None, comps=None): """De-project from a map, returning a Signal-like object. From f0ef799cac3e414eeee6fb6d7935204d6ea56466 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Thu, 16 Jan 2025 09:58:02 -0800 Subject: [PATCH 12/13] Fixed last-minute typo --- python/proj/coords.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/proj/coords.py b/python/proj/coords.py index 1fcb2ae1..dbcf4985 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -231,7 +231,7 @@ def coords(self, fplane=None, output=None): fplane = FocalPlane.boresight() if output is not None: output = output[None] - output = p.coords(self.Q, fplane.quats), output) + output = p.coords(self.Q, fplane.quats, output) if collapse: output = output[0] return output From 263d37821d17c3413a0ada9753bc7c9dfc654024 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sigurd=20Kirkevold=20N=C3=A6ss?= Date: Fri, 17 Jan 2025 02:04:37 -0800 Subject: [PATCH 13/13] Addressed last round of commments --- python/proj/coords.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/python/proj/coords.py b/python/proj/coords.py index dbcf4985..11d8bf95 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -240,11 +240,11 @@ class FocalPlane: """This class represents the detector positions and intensity and polarization responses in the focal plane. - Members: + Attributes: quats: G3VectorQuat representing the pointing quaternions for each detector. Can be turned into a numpy array of coefficients with np.array(). Or call .coeffs() to get them directly. - resps: Array of float64 with shape [ndet,2] representing the + resps: Array of float32 with shape [ndet,2] representing the total intensity and polarization responsivity of each detector ndet: The number of detectors (read-only) @@ -300,9 +300,13 @@ def boresight(cls): @classmethod def from_xieta(cls, *args, **kwargs): """ - def from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False): + ``from_xieta(cls, xi, eta, gamma=0, T=1, P=1, Q=1, U=0, hwp=False)`` Construct a FocalPlane from focal plane coordinates (xi,eta). + For backwards compatibility, the signature + ``from_xieta(names, xi, eta, gamma=0)`` is also supported; but + this will be removed in the future. + These are Cartesian projection plane coordinates. When looking at the sky along the telescope boresight, xi is parallel to increasing azimuth and eta is parallel to increasing elevation.