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/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 4ddb42e2..9401ada6 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 @@ -110,26 +110,57 @@ 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:: - - >>> fp['c'] +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] 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. 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:: + + 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]) + 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 * fp.quats[0] # Decompose q_total into lon, lat, roll angles ra, dec, gamma = so3g.proj.quat.decompose_lonlat(q_total) @@ -137,21 +168,19 @@ 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 FocalPlane object as first argument:: >>> csl.coords(fp) - 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, + [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]])] + +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))``. @@ -215,17 +244,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`:: @@ -240,9 +269,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 @@ -251,7 +280,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 @@ -264,8 +293,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 @@ -289,7 +318,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]]) @@ -352,7 +381,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/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 7ffb5fd6..11d8bf95 100644 --- a/python/proj/coords.py +++ b/python/proj/coords.py @@ -208,93 +208,220 @@ 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. - output: An optional structure for holding the results. For + 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 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.boresight() 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, 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. + 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 float32 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 + # for this when we use FocalPlane later + 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 + 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") + # 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""" + return cls(quats=[[1,0,0,0]]) + # FIXME: old sotodlib compat - expand to actual argument list later @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. + def from_xieta(cls, *args, **kwargs): + """ + ``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). - 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. + For backwards compatibility, the signature + ``from_xieta(names, xi, eta, gamma=0)`` is also supported; but + this will be removed in the future. - """ - qs = quat.rotation_xieta(np.asarray(xi), np.asarray(eta), np.asarray(gamma)) - return cls([(n,q) for n, q in zip(names, qs)]) + 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. 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``. + """ + # 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(*args, **kwargs) + gamma = gamma + np.arctan2(U,Q)/2 + P = P * (Q**2+U**2)**0.5 + 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 + # 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.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 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 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 + 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(self.quats, self.resps): + yield q, resp + # FIXME: old sotodlib compat - remove later + @staticmethod + def _xieta_compat(*args, **kwargs): + # Accept the alternative format (names,xi,eta,gamma) + 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 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 + if name in self._detmap: + i = self._detmap[i] + 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.resps = np.concatenate([self.resps,np.ones((1,2),np.float32)]) class Assembly: """This class groups together boresight and detector offset @@ -305,12 +432,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,34 +444,26 @@ 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) + self = cls() 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 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.dets = [np.array([1., 0, 0, 0])] + self.fplane = FocalPlane.boresight() return self + # FIXME: old sotodlib compat - remove later + @property + def dets(self): return self.fplane 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 3b355ccd..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.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 @@ -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.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. @@ -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.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. @@ -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.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,) + @@ -272,8 +272,8 @@ 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.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, @@ -298,8 +298,8 @@ 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.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): @@ -318,8 +318,8 @@ 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.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): @@ -362,11 +362,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'] @@ -374,9 +374,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) @@ -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.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): @@ -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.dets)) + 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.dets, 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.dets, 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. diff --git a/src/Projection.cxx b/src/Projection.cxx index 4ab1789d..96fc7a77 100644 --- a/src/Projection.cxx +++ b/src/Projection.cxx @@ -22,6 +22,7 @@ using namespace std; #include "exceptions.h" #include "so_linterp.h" +#include #include "healpix_bare.c" // TRIG_TABLE_SIZE @@ -122,6 +123,31 @@ 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 +// ----------------------------------- +// 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. +// +// 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. +// +// 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. +// +// Finally, this Response struct is passed to spin_proj_factors, e.g. +// +// 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; class ProjFlat; @@ -1094,35 +1120,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); } @@ -1307,7 +1345,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(); @@ -1320,9 +1358,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) { @@ -1332,12 +1371,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] + @@ -1642,7 +1682,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(); @@ -1658,6 +1698,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) { @@ -1666,10 +1708,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) @@ -1684,6 +1727,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, @@ -1701,13 +1745,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. @@ -1723,6 +1768,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) @@ -1739,12 +1785,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? @@ -1831,8 +1878,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(); @@ -1852,6 +1899,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 @@ -1866,15 +1915,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(); @@ -1894,6 +1943,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 @@ -1908,7 +1959,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; } 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..17bd261c 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,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.dets[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 a7eb88c0..b407ade8 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) @@ -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.dets[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