diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3683b6ad..e0fc9f25 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -56,6 +56,10 @@ jobs: pip install --no-build-isolation --no-deps -e . cd .. + - name: Try newer healpy + run: | + pip install healpy>=1.17.3 --upgrade + - name: Cache data id: cache-data uses: actions/cache@v4 diff --git a/imsim/treerings.py b/imsim/treerings.py index a10939cf..fb822aac 100644 --- a/imsim/treerings.py +++ b/imsim/treerings.py @@ -6,10 +6,49 @@ from .meta_data import data_dir import warnings + class TreeRingsError(Exception): pass -class TreeRings(): + +class TreeRingRadialFunction: + """ + Radial Function describing tree rings in a CCD. + """ + def __init__(self, info_block): + """ + Parameters + ---------- + info_block : list of text lines + This is the block of data describing the radial tree ring profile for a + specific CCD. + """ + items = info_block[1].split() + self.A = float(items[6]) + self.B = float(items[7]) + self.cfreqs, self.cphases, self.sfreqs, self.sphases = np.genfromtxt(info_block[3:]).T + + def __call__(self, r): + """ + This function defines the tree ring distortion of the pixels as + a radial function. + + Parameters + ---------- + r: float + Radial coordinate from the center of the tree ring structure + in units of pixels. + """ + centroid_shift = 0.0 + for j, fval in enumerate(self.cfreqs): + centroid_shift += np.sin(2*np.pi*(r/fval)+self.cphases[j]) * fval / (2.0*np.pi) + for j, fval in enumerate(self.sfreqs): + centroid_shift += -np.cos(2*np.pi*(r/fval)+self.sphases[j]) * fval / (2.0*np.pi) + centroid_shift *= (self.A + self.B * r**4) * .01 # 0.01 factor is because data is in percent + return centroid_shift + + +class TreeRings: """ # Craig Lage UC Davis 16-Mar-18; cslage@ucdavis.edu # This function returns a tree ring model drawn from an analytical function that was @@ -20,7 +59,8 @@ class TreeRings(): # amplitude 10X greater that the 60% of the sensors that have 'good' tree rings. """ _req_params = {'file_name' : str} - _opt_params = {'only_dets' : list} + _opt_params = {'only_dets' : list, + 'defer_load' : bool} def __init__(self, file_name, only_dets=None, logger=None, defer_load=True): """ @@ -38,15 +78,18 @@ def __init__(self, file_name, only_dets=None, logger=None, defer_load=True): raise OSError("TreeRing file %s not found"%file_name) self.only_dets = only_dets self.numfreqs = 20 # Number of spatial frequencies - self.cfreqs = np.zeros([self.numfreqs]) # Cosine frequencies - self.cphases = np.zeros([self.numfreqs]) - self.sfreqs = np.zeros([self.numfreqs]) # Sine frequencies - self.sphases = np.zeros([self.numfreqs]) self.r_max = 8000.0 # Maximum extent of tree ring function in pixels dr = 3.0 # Step size of tree ring function in pixels self.npoints = int(self.r_max / dr) + 1 # Number of points in tree ring function logger.warning("TreeRing file %s will be used.", self.file_name) + self._read_info_blocks() + # Check for entries in only_dets that are not in info_blocks: + if only_dets: + missing_dets = set(only_dets).difference(self.info_blocks) + if missing_dets: + logger.info("Requested det_names that are not in the " + "tree ring info file: %s", missing_dets) # Make a dict indexed by det_name (a string that looks like Rxx_Syy) self.info = {} if not defer_load: @@ -54,36 +97,48 @@ def __init__(self, file_name, only_dets=None, logger=None, defer_load=True): logger.info("Reading in det_names: %s", only_dets) self.fill_dict(only_dets=only_dets) - def fill_dict(self, only_dets=None): - """Read the tree rings file and save the information therein to a dict, self.info. + def _read_info_blocks(self): + """ + Read the tree ring data for all of the CCDs in the info file, filling a + dictionary keyed by detector name. """ - with open(self.file_name, 'r') as input_: - lines = input_.readlines() + with open(self.file_name, 'r') as fobj: + lines = fobj.readlines() + block_size = self.numfreqs + 3 + self.info_blocks = {} + num_blocks = len(lines) // block_size + for iblock in range(num_blocks): + imin = iblock*block_size + imax = imin + block_size + block = lines[imin:imax] + items = block[1].split() + det_name = "R%s%s_S%s%s" % (tuple(items[:4])) + self.info_blocks[det_name] = block + def fill_dict(self, only_dets=None): + """ + Fill the self.info dictionary with the tree ring model for the detectors in + the info file, restricting to those listed in only_dets, if provided. + """ xCenterPix = 2048.5 yCenterPix = 2048.5 - for i, line in enumerate(lines): - if not line.startswith('Rx'): continue - items = lines[i+1].split() + if only_dets is None: + # Process all detectors in the tree ring info file. + only_dets = self.info_blocks.keys() - det_name = "R%s%s_S%s%s"%(tuple(items[:4])) - if only_dets and det_name not in only_dets: continue + for det_name in only_dets: + if det_name not in self.info_blocks: + continue + info_block = self.info_blocks[det_name] + items = info_block[1].split() Cx = float(items[4]) + xCenterPix Cy = float(items[5]) + yCenterPix center = galsim.PositionD(Cx, Cy) - # Temporary instance variables used by tree_ring_radial_function - self.A = float(items[6]) - self.B = float(items[7]) - for j in range(self.numfreqs): - freqitems = lines[i + 3 + j].split() - self.cfreqs[j] = float(freqitems[0]) - self.cphases[j] = float(freqitems[1]) - self.sfreqs[j] = float(freqitems[2]) - self.sphases[j] = float(freqitems[3]) - func = galsim.LookupTable.from_func(self.tree_ring_radial_function, + tree_ring_radial_function = TreeRingRadialFunction(info_block) + func = galsim.LookupTable.from_func(tree_ring_radial_function, x_min=0.0, x_max=self.r_max, npoints=self.npoints) self.info[det_name] = (center, func) @@ -106,25 +161,6 @@ def get_func(self, det_name): warnings.warn("No treering information available for %s. Setting treering_func to None." % det_name) return None - def tree_ring_radial_function(self, r): - """ - This function defines the tree ring distortion of the pixels as - a radial function. - - Parameters - ---------- - r: float - Radial coordinate from the center of the tree ring structure - in units of pixels. - """ - centroid_shift = 0.0 - for j, fval in enumerate(self.cfreqs): - centroid_shift += np.sin(2*np.pi*(r/fval)+self.cphases[j]) * fval / (2.0*np.pi) - for j, fval in enumerate(self.sfreqs): - centroid_shift += -np.cos(2*np.pi*(r/fval)+self.sphases[j]) * fval / (2.0*np.pi) - centroid_shift *= (self.A + self.B * r**4) * .01 # 0.01 factor is because data is in percent - return centroid_shift - def TreeRingCenter(config, base, value_type): """Return the tree ring center for the current det_name. """ diff --git a/tests/test_tree_rings.py b/tests/test_tree_rings.py index a99798c7..36b8bb66 100644 --- a/tests/test_tree_rings.py +++ b/tests/test_tree_rings.py @@ -73,6 +73,16 @@ def test_deferred_read(self): self.assertEqual(len(tree_rings.info), 2) self.assertIn(det_name, tree_rings.info) + def test_unhandled_det_names(self): + """Test that unrecognized detectors in only_dets are handled + without raising errors.""" + # There is no detector named R44_S12 in LSSTCam. + only_dets = ['R22_S12', 'R44_S12'] + tree_rings = imsim.TreeRings(self.tr_filename, only_dets=only_dets) + for det_name in only_dets: + _ = tree_rings.get_center(det_name) + _ = tree_rings.get_func(det_name) + if __name__ == '__main__': unittest.main()