22932166231

Committed 11 Mar 2026 01:25AM UTC coverage: 56.018%. Remained the same

Build # 22932166231

Build Type

push

github

Committed by

mhasself

Commit Message

Fix domdir thread assignement when tiled and all tiles empty.

Run Details

7 of 7 new or added lines in 1 file covered. (100.0%)

29 existing lines in 1 file now uncovered.

1345 of 2401 relevant lines covered (56.02%)

0.56 hits per line

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

65.98

/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
from .ranges import RangesMatrix
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.
    if len(pmat.active_tiles) > 0:
        T, Q, U = np.sum([_m.reshape((3, -1)).sum(axis=-1)
                          for _m in scan_maps if _m is not None], axis=0)
    else:
        T, Q, U = 0., 0., 0.

    if T == 0:
        # No hits -- plan to do nothing, but use n_threads to avoid
        # needing a special test case.
        n_dets, n_samps = len(asm_full.dets), len(asm_full.Q)
        return [RangesMatrix.zeros(shape=(n_threads, n_dets, n_samps))]

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

simonsobs / so3g / 22932166231

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