12832805671

Committed 17 Jan 2025 04:16PM UTC coverage: 50.45% (+0.4%) from 50.055%

Build # 12832805671

Build Type

push

github

Committed by

web-flow

Commit Message

Stokes response (#193)

* Plans for arbitrary detector response stuff

* Detector response now supported in c++, and represented in a redone FocalPlane class in python

* Arbitrary response support

* Passes tests after adding .items(), updating xieta arguments and replacing asm.dets with asm.fplane.quats

* Temporary compatibility with old sotodlib

* Default FocalPlane is not empty, so make a separate constructor for making one with just the boresight

* Better docstrings and comments

* Working on new documentation

* More documentation. Updated numbers in examples which seem to have rotted at some point

* Matthews comments, plus go to G3VectorQuat instead of numpy array

* Ready for new review

* Fixed last-minute typo

* Addressed last round of commments

Run Details

69 of 87 new or added lines in 3 files covered. (79.31%)

1 existing line in 1 file now uncovered.

1401 of 2777 relevant lines covered (50.45%)

0.5 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

63.74

/python/proj/mapthreads.py

"""This submodule is for functions that produce "thread assignments",
i.e. disjoint RangesMatrix objects with shape (n_threads, n_dets,
n_samps) that are suitable for the threads= argument in Projection
code that projects from time to map space using OpenMP
parallelization.

"""

import so3g
import numpy as np

def get_num_threads(n_threads=None):
    """Utility function for computing n_threads.  If n_threads is not
    None, it is returned directly.  But if it is None, then the OpenMP
    thread count is returned.  Uses so3g.useful_info().

    """
    if n_threads is None:
        return so3g.useful_info()['omp_num_threads']
    return n_threads

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.

    The dominant scan direction is first determined, which requires
    creating a 3-component map.  This projection operation can't be
    parallelized safely so it is wise to pass in a decimated detector
    set through offs_rep.  The second step is to assign pixels to
    threads; to do this a signal-sized array is projected from map
    space.  In the present implementation this uses all the detectors
    (but takes advantage of threads).  In principle this step could be
    done with a decimated detector set, though with care that the
    decimated subset covers the array well, spatially.

    Arguments:
      sight (CelestialSightLine): The boresight 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
        tiles.
      active_tiles (list): The list of active tiles.  If None, this
        will be computed along the way.
      n_threads (int): The number of threads to target (defaults to
        OMP_NUM_THREADS).
      fplane_rep (FocalPlane): A representative set of
        detectors, for determining scan direction (if not present then
        fplane is used for this).
      plot_prefix (str): If present, pylab will be imported and some
        plots saved with this name as prefix.

    Returns:
      Ranges matrix with shape (n_thread, n_det, n_samp), that can be
      passes as the "threads" argument to Projectionist routines.

    """
    if plot_prefix:
        import pylab as pl
        if tile_shape is None:
            tile_iter = lambda maps: [('', maps[0])]
        else:
            tile_iter = lambda maps: [('_tile%04i' % i, _m)
                                      for i, _m in enumerate(maps)
                                      if _m is not None]

    n_threads = get_num_threads(n_threads)
    if fplane_rep is None:
        fplane_rep = fplane

    if tile_shape is None:
        # Let's pretend it is, though; this simplifies logic below and
        # doesn't cost much.
        tile_shape = shape
        active_tiles = [0]

    # The full assembly, for later.
    asm_full = so3g.proj.Assembly.attach(sight, fplane)

    # Get a Projectionist -- note it can be used with full or
    # representative assembly.
    pmat = so3g.proj.wcs.Projectionist.for_tiled(
        shape, wcs, tile_shape=tile_shape, active_tiles=active_tiles
    )
    if active_tiles is None:
        # This is OMPed, but it's still better to pass in active_tiles
        # if you've already computed it somehow.
        pmat.active_tiles = pmat.get_active_tiles(asm_full, assign=True)['active_tiles']
    tiling = pmat.tiling

    # 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(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
    # that +Q is parallel to map columns, and U increases along the
    # map diagonal.
    T, Q, U = np.sum([_m.reshape((3, -1)).sum(axis=-1)
                      for _m in scan_maps if _m is not None], axis=0)
    phi = np.arctan2(U, Q) / 2

    if plot_prefix:
        text = 'Qf=%.2f Uf=%.2f phi=%.1f deg' % (Q/T, U/T, phi / so3g.proj.DEG)
        for label, _m in tile_iter(scan_maps):
            for i in range(3):
                pl.imshow(_m[i], origin='lower')
                pl.title(label + ' '  + 'TQU'[i] + ' ' + text)
                pl.colorbar()
                pl.savefig(plot_prefix + '%02i%s%s.png' % (i, 'TQU'[i], label))
                pl.clf()

    # Now paint a map with a single parameter such that contours are
    # parallel to the scan direction.
    idx_maps = pmat.zeros((1,))
    lims = None
    for t in pmat.active_tiles:
        y0, x0 = tiling.tile_offset(t)
        ny, nx = tiling.tile_dims(t)
        y, x = y0 + np.arange(ny), x0 + np.arange(nx)
        idx_maps[t] = y[:,None]*np.sin(phi) - x[None,:]*np.cos(phi)
        _lo, _hi = idx_maps[t].min(), idx_maps[t].max()
        if lims is None:
            lims = _lo, _hi
        else:
            lims = min(_lo, lims[0]), max(_hi, lims[1])

    if plot_prefix:
        for label, _m in tile_iter(idx_maps):
            pl.imshow(_m, origin='lower')
            pl.colorbar()
            pl.savefig(plot_prefix + '10param%s.png' % label)
            pl.clf()

    # Histogram the parameter, weighted by the weights map -- this
    # tells us the number of hits for ranges of the parameter value.
    n_bins = 200
    bins = np.linspace(lims[0], lims[1], n_bins+1)
    H = np.zeros(n_bins)
    for t in pmat.active_tiles:
        H += np.histogram(idx_maps[t].ravel(), bins=bins,
                          weights=scan_maps[t][0].ravel())[0]
    del scan_maps

    # Turn histogram into a CDF.  The CDF is non-decreasing, but could
    # have some flat parts in it, so bias the whole curve to guarantee
    # it's strictly increasing.
    H = np.hstack((0, H.cumsum()))
    bias = .00001
    H = H / H.max() * (1-bias) + np.linspace(0, bias, len(H))

    # Create n_threads "super_bins" that contain roughly equal weight.
    xs = np.linspace(0, 1, n_threads + 1)
    ys = np.interp(xs, H, bins)
    superbins = np.hstack((bins[0], ys[1:-1], bins[-1]))

    # Create maps where the pixel value is a superbin index.
    for t in pmat.active_tiles:
        temp = np.zeros(idx_maps[t].shape)
        for i in range(n_threads):
            s = (superbins[i] <= idx_maps[t])*(idx_maps[t] < superbins[i+1])
            temp[s] = i
        idx_maps[t][:] = temp
        idx_maps[t].shape = (1, ) + tuple(idx_maps[t].shape)

    # Turn the superbin index maps into thread assignments.
    threads = pmat.assign_threads_from_map(asm_full, idx_maps, n_threads=n_threads)

    if plot_prefix:
        pl.plot(bins, H)
        for _x, _y in zip(superbins[1:-1], xs[1:-1]):
            pl.plot([_x, _x], [-1, _y], c='k')
            pl.plot([superbins[0], superbins[-1]], [_y, _y], c='gray')
        pl.xlim(superbins[0], superbins[-1])
        pl.ylim(-.01, 1.01)
        pl.savefig(plot_prefix + '20histo.png')
        pl.clf()

        for label, _m in tile_iter(idx_maps):
            pl.imshow(_m[0], origin='lower')
            pl.colorbar()
            pl.savefig(plot_prefix + '30thread%s.png' % label)
            pl.clf()

    return threads


1	"""This submodule is for functions that produce "thread assignments",
2	i.e. disjoint RangesMatrix objects with shape (n_threads, n_dets,
3	n_samps) that are suitable for the threads= argument in Projection
4	code that projects from time to map space using OpenMP
5	parallelization.
6
7	"""
8
9	import so3g	1✔
10	import numpy as np	1✔
11
12	def get_num_threads(n_threads=None):	1✔
13	"""Utility function for computing n_threads. If n_threads is not
14	None, it is returned directly. But if it is None, then the OpenMP
15	thread count is returned. Uses so3g.useful_info().
16
17	"""
18	if n_threads is None:	1✔
19	return so3g.useful_info()['omp_num_threads']	×
20	return n_threads	1✔
21
22	def get_threads_domdir(sight, fplane, shape, wcs, tile_shape=None,	1✔
23	active_tiles=None, n_threads=None, fplane_rep=None,
24	plot_prefix=None):
25	"""Assign samples to threads according to the dominant scan
26	direction.
27
28	The dominant scan direction is first determined, which requires
29	creating a 3-component map. This projection operation can't be
30	parallelized safely so it is wise to pass in a decimated detector
31	set through offs_rep. The second step is to assign pixels to
32	threads; to do this a signal-sized array is projected from map
33	space. In the present implementation this uses all the detectors
34	(but takes advantage of threads). In principle this step could be
35	done with a decimated detector set, though with care that the
36	decimated subset covers the array well, spatially.
37
38	Arguments:
39	sight (CelestialSightLine): The boresight pointing.
40	fplane (FocalPlane): The detector pointing
41	shape (tuple): The map shape.
42	wcs (wcs): The map WCS.
43	tile_shape (tuple): The tile shape, if this should be done using
44	tiles.
45	active_tiles (list): The list of active tiles. If None, this
46	will be computed along the way.
47	n_threads (int): The number of threads to target (defaults to
48	OMP_NUM_THREADS).
49	fplane_rep (FocalPlane): A representative set of
50	detectors, for determining scan direction (if not present then
51	fplane is used for this).
52	plot_prefix (str): If present, pylab will be imported and some
53	plots saved with this name as prefix.
54
55	Returns:
56	Ranges matrix with shape (n_thread, n_det, n_samp), that can be
57	passes as the "threads" argument to Projectionist routines.
58
59	"""
60	if plot_prefix:	1✔
61	import pylab as pl	×
62	if tile_shape is None:	×
63	tile_iter = lambda maps: [('', maps[0])]	×
64	else:
65	tile_iter = lambda maps: [('_tile%04i' % i, _m)	×
66	for i, _m in enumerate(maps)
67	if _m is not None]
68
69	n_threads = get_num_threads(n_threads)	1✔
70	if fplane_rep is None:	1✔
NEW 71	fplane_rep = fplane	×
72
73	if tile_shape is None:	1✔
74	# Let's pretend it is, though; this simplifies logic below and
75	# doesn't cost much.
76	tile_shape = shape	1✔
77	active_tiles = [0]	1✔
78
79	# The full assembly, for later.
80	asm_full = so3g.proj.Assembly.attach(sight, fplane)	1✔
81
82	# Get a Projectionist -- note it can be used with full or
83	# representative assembly.
84	pmat = so3g.proj.wcs.Projectionist.for_tiled(	1✔
85	shape, wcs, tile_shape=tile_shape, active_tiles=active_tiles
86	)
87	if active_tiles is None:	1✔
88	# This is OMPed, but it's still better to pass in active_tiles
89	# if you've already computed it somehow.
90	pmat.active_tiles = pmat.get_active_tiles(asm_full, assign=True)['active_tiles']	×
91	tiling = pmat.tiling	1✔
92
93	# For the scan direction map, use the "representative" subset
94	# detectors, with polarization direction aligned parallel to
95	# elevation.
96	xi, eta, gamma = so3g.proj.quat.decompose_xieta(fplane_rep.quats)	1✔
97	fplane_xl = so3g.proj.FocalPlane.from_xieta(xi, eta, gamma0+90so3g.proj.DEG)	1✔
98	asm_rep = so3g.proj.Assembly.attach(sight, fplane_xl)	1✔
99	sig = np.ones((fplane_xl.ndet, len(asm_rep.Q)), dtype='float32')	1✔
100	scan_maps = pmat.to_map(sig, asm_rep, comps='TQU')	1✔
101
102	# Compute the scan angle, based on Q and U weights. This assumes
103	# that +Q is parallel to map columns, and U increases along the
104	# map diagonal.
105	T, Q, U = np.sum([_m.reshape((3, -1)).sum(axis=-1)	1✔
106	for _m in scan_maps if _m is not None], axis=0)
107	phi = np.arctan2(U, Q) / 2	1✔
108
109	if plot_prefix:	1✔
110	text = 'Qf=%.2f Uf=%.2f phi=%.1f deg' % (Q/T, U/T, phi / so3g.proj.DEG)	×
111	for label, _m in tile_iter(scan_maps):	×
112	for i in range(3):	×
113	pl.imshow(_m[i], origin='lower')	×
114	pl.title(label + ' ' + 'TQU'[i] + ' ' + text)	×
115	pl.colorbar()	×
116	pl.savefig(plot_prefix + '%02i%s%s.png' % (i, 'TQU'[i], label))	×
117	pl.clf()	×
118
119	# Now paint a map with a single parameter such that contours are
120	# parallel to the scan direction.
121	idx_maps = pmat.zeros((1,))	1✔
122	lims = None	1✔
123	for t in pmat.active_tiles:	1✔
124	y0, x0 = tiling.tile_offset(t)	1✔
125	ny, nx = tiling.tile_dims(t)	1✔
126	y, x = y0 + np.arange(ny), x0 + np.arange(nx)	1✔
127	idx_maps[t] = y[:,None]np.sin(phi) - x[None,:]np.cos(phi)	1✔
128	_lo, _hi = idx_maps[t].min(), idx_maps[t].max()	1✔
129	if lims is None:	1✔
130	lims = _lo, _hi	1✔
131	else:
132	lims = min(_lo, lims[0]), max(_hi, lims[1])	1✔
133
134	if plot_prefix:	1✔
135	for label, _m in tile_iter(idx_maps):	×
136	pl.imshow(_m, origin='lower')	×
137	pl.colorbar()	×
138	pl.savefig(plot_prefix + '10param%s.png' % label)	×
139	pl.clf()	×
140
141	# Histogram the parameter, weighted by the weights map -- this
142	# tells us the number of hits for ranges of the parameter value.
143	n_bins = 200	1✔
144	bins = np.linspace(lims[0], lims[1], n_bins+1)	1✔
145	H = np.zeros(n_bins)	1✔
146	for t in pmat.active_tiles:	1✔
147	H += np.histogram(idx_maps[t].ravel(), bins=bins,	1✔
148	weights=scan_maps[t][0].ravel())[0]
149	del scan_maps	1✔
150
151	# Turn histogram into a CDF. The CDF is non-decreasing, but could
152	# have some flat parts in it, so bias the whole curve to guarantee
153	# it's strictly increasing.
154	H = np.hstack((0, H.cumsum()))	1✔
155	bias = .00001	1✔
156	H = H / H.max() * (1-bias) + np.linspace(0, bias, len(H))	1✔
157
158	# Create n_threads "super_bins" that contain roughly equal weight.
159	xs = np.linspace(0, 1, n_threads + 1)	1✔
160	ys = np.interp(xs, H, bins)	1✔
161	superbins = np.hstack((bins[0], ys[1:-1], bins[-1]))	1✔
162
163	# Create maps where the pixel value is a superbin index.
164	for t in pmat.active_tiles:	1✔
165	temp = np.zeros(idx_maps[t].shape)	1✔
166	for i in range(n_threads):	1✔
167	s = (superbins[i] <= idx_maps[t])*(idx_maps[t] < superbins[i+1])	1✔
168	temp[s] = i	1✔
169	idx_maps[t][:] = temp	1✔
170	idx_maps[t].shape = (1, ) + tuple(idx_maps[t].shape)	1✔
171
172	# Turn the superbin index maps into thread assignments.
173	threads = pmat.assign_threads_from_map(asm_full, idx_maps, n_threads=n_threads)	1✔
174
175	if plot_prefix:	1✔
176	pl.plot(bins, H)	×
177	for _x, _y in zip(superbins[1:-1], xs[1:-1]):	×
178	pl.plot([_x, _x], [-1, _y], c='k')	×
179	pl.plot([superbins[0], superbins[-1]], [_y, _y], c='gray')	×
180	pl.xlim(superbins[0], superbins[-1])	×
181	pl.ylim(-.01, 1.01)	×
182	pl.savefig(plot_prefix + '20histo.png')	×
183	pl.clf()	×
184
185	for label, _m in tile_iter(idx_maps):	×
186	pl.imshow(_m[0], origin='lower')	×
187	pl.colorbar()	×
188	pl.savefig(plot_prefix + '30thread%s.png' % label)	×
189	pl.clf()	×
190
191	return threads	1✔
192

simonsobs / so3g / 12832805671

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous