diff --git a/py/fiberassign/test/simulate.py b/py/fiberassign/test/simulate.py index 236e0b94..c7ba3a26 100644 --- a/py/fiberassign/test/simulate.py +++ b/py/fiberassign/test/simulate.py @@ -21,7 +21,8 @@ FIBER_STATE_BROKEN) from fiberassign.targets import (TARGET_TYPE_SCIENCE, TARGET_TYPE_SKY, - TARGET_TYPE_STANDARD) + TARGET_TYPE_SUPPSKY, + TARGET_TYPE_SUPPSKY, TARGET_TYPE_STANDARD) def sim_data_dir(): @@ -167,6 +168,7 @@ def sim_targets(path, tgtype, tgoffset, density=5000.0): fdata["SUBPRIORITY"] = np.random.uniform(low=0.0, high=1.0, size=ntarget) sky_mask = desi_mask["SKY"].mask + suppsky_mask = desi_mask["SUPP_SKY"].mask std_mask = desi_mask["STD_BRIGHT"].mask # We could be fancier and set the DESI_TARGET bits for science targets # to exactly match the priority class. For now we just set a single @@ -180,6 +182,9 @@ def sim_targets(path, tgtype, tgoffset, density=5000.0): elif tgtype == TARGET_TYPE_STANDARD: fdata["PRIORITY"] = 1500 * np.ones(ntarget, dtype=np.int32) fdata["DESI_TARGET"] |= std_mask + elif tgtype == TARGET_TYPE_SUPPSKY: + fdata["PRIORITY"] = np.zeros(ntarget, dtype=np.int32) + fdata["DESI_TARGET"] |= suppsky_mask else: fdata["DESI_TARGET"] |= sci_mask # These are the fractions of each target type: diff --git a/py/fiberassign/test/test_assign.py b/py/fiberassign/test/test_assign.py index 91f69648..a51138cf 100644 --- a/py/fiberassign/test/test_assign.py +++ b/py/fiberassign/test/test_assign.py @@ -24,6 +24,7 @@ from fiberassign.tiles import load_tiles, Tiles from fiberassign.targets import (TARGET_TYPE_SCIENCE, TARGET_TYPE_SKY, + TARGET_TYPE_SUPPSKY, TARGET_TYPE_STANDARD, TARGET_TYPE_SAFE, Targets, TargetsAvailable, TargetTree, LocationsAvailable, load_target_file) @@ -52,6 +53,10 @@ class TestAssign(unittest.TestCase): def setUp(self): + self.density_science = 5000 + self.density_standards = 5000 + self.density_sky = 10 + self.density_suppsky = 5000 pass def tearDown(self): @@ -62,14 +67,41 @@ def test_io(self): input_mtl = os.path.join(test_dir, "mtl.fits") input_std = os.path.join(test_dir, "standards.fits") input_sky = os.path.join(test_dir, "sky.fits") - nscience = sim_targets(input_mtl, TARGET_TYPE_SCIENCE, 0) - nstd = sim_targets(input_std, TARGET_TYPE_STANDARD, nscience) - nsky = sim_targets(input_sky, TARGET_TYPE_SKY, (nscience + nstd)) + input_suppsky = os.path.join(test_dir, "suppsky.fits") + tgoff = 0 + nscience = sim_targets( + input_mtl, + TARGET_TYPE_SCIENCE, + tgoff, + density=self.density_science + ) + tgoff += nscience + nstd = sim_targets( + input_std, + TARGET_TYPE_STANDARD, + tgoff, + density=self.density_standards + ) + tgoff += nstd + nsky = sim_targets( + input_sky, + TARGET_TYPE_SKY, + tgoff, + density=self.density_sky + ) + tgoff += nsky + nsuppsky = sim_targets( + input_suppsky, + TARGET_TYPE_SUPPSKY, + tgoff, + density=self.density_suppsky + ) tgs = Targets() load_target_file(tgs, input_mtl) load_target_file(tgs, input_std) load_target_file(tgs, input_sky) + load_target_file(tgs, input_suppsky) # Create a hierarchical triangle mesh lookup of the targets positions tree = TargetTree(tgs, 0.01) @@ -103,16 +135,19 @@ def test_io(self): write_assignment_fits(tiles, asgn, out_dir=test_dir, out_prefix="full_", all_targets=True) + plotpetals = [0] + # plotpetals = None + plot_tiles(hw, tiles, result_dir=test_dir, result_prefix="basic_", plot_dir=test_dir, plot_prefix="basic_", - result_split_dir=False, petals=[0], + result_split_dir=False, petals=plotpetals, serial=True) plot_tiles(hw, tiles, result_dir=test_dir, result_prefix="full_", plot_dir=test_dir, plot_prefix="full_", - result_split_dir=False, petals=[0], + result_split_dir=False, petals=plotpetals, serial=True) target_files = [ @@ -249,13 +284,13 @@ def test_io(self): plot_tiles(hw, tiles, result_dir=test_dir, result_prefix="basic_tile-", plot_dir=test_dir, plot_prefix="basic_tile-", - result_split_dir=False, petals=[0], + result_split_dir=False, petals=plotpetals, serial=True) plot_tiles(hw, tiles, result_dir=test_dir, result_prefix="full_tile-", plot_dir=test_dir, plot_prefix="full_tile-", - result_split_dir=False, petals=[0], + result_split_dir=False, petals=plotpetals, serial=True) return @@ -265,14 +300,41 @@ def test_full(self): input_mtl = os.path.join(test_dir, "mtl.fits") input_std = os.path.join(test_dir, "standards.fits") input_sky = os.path.join(test_dir, "sky.fits") - nscience = sim_targets(input_mtl, TARGET_TYPE_SCIENCE, 0) - nstd = sim_targets(input_std, TARGET_TYPE_STANDARD, nscience) - nsky = sim_targets(input_sky, TARGET_TYPE_SKY, (nscience + nstd)) + input_suppsky = os.path.join(test_dir, "suppsky.fits") + tgoff = 0 + nscience = sim_targets( + input_mtl, + TARGET_TYPE_SCIENCE, + tgoff, + density=self.density_science + ) + tgoff += nscience + nstd = sim_targets( + input_std, + TARGET_TYPE_STANDARD, + tgoff, + density=self.density_standards + ) + tgoff += nstd + nsky = sim_targets( + input_sky, + TARGET_TYPE_SKY, + tgoff, + density=self.density_sky + ) + tgoff += nsky + nsuppsky = sim_targets( + input_suppsky, + TARGET_TYPE_SUPPSKY, + tgoff, + density=self.density_suppsky + ) tgs = Targets() load_target_file(tgs, input_mtl) load_target_file(tgs, input_std) load_target_file(tgs, input_sky) + load_target_file(tgs, input_suppsky) # Create a hierarchical triangle mesh lookup of the targets positions tree = TargetTree(tgs, 0.01) @@ -310,14 +372,21 @@ def test_full(self): asgn.assign_unused(TARGET_TYPE_SKY, 40) asgn.assign_force(TARGET_TYPE_SKY, 40) + # Use supplemental sky to meet our requirements + asgn.assign_unused(TARGET_TYPE_SUPPSKY, 40) + asgn.assign_force(TARGET_TYPE_SUPPSKY, 40) + # If there are any unassigned fibers, try to place them somewhere. asgn.assign_unused(TARGET_TYPE_SCIENCE) asgn.assign_unused(TARGET_TYPE_SKY) + asgn.assign_unused(TARGET_TYPE_SUPPSKY) write_assignment_fits(tiles, asgn, out_dir=test_dir, all_targets=True) + # plotpetals = [0] + plotpetals = None plot_tiles(hw, tiles, result_dir=test_dir, plot_dir=test_dir, - real_shapes=True, petals=[4, 6], serial=True) + real_shapes=True, petals=plotpetals, serial=True) qa_tiles(hw, tiles, result_dir=test_dir) @@ -328,7 +397,9 @@ def test_full(self): for tile, props in qadata.items(): self.assertEqual(4495, props["assign_science"]) self.assertEqual(100, props["assign_std"]) - self.assertEqual(400, props["assign_sky"]) + self.assertEqual(400, ( + props["assign_sky"] + props["assign_suppsky"] + )) plot_qa(qadata, os.path.join(test_dir, "qa"), outformat="pdf", labels=True) @@ -340,16 +411,41 @@ def test_cli(self): input_mtl = os.path.join(test_dir, "mtl.fits") input_std = os.path.join(test_dir, "standards.fits") input_sky = os.path.join(test_dir, "sky.fits") - - nscience = sim_targets(input_mtl, TARGET_TYPE_SCIENCE, 0) - nstd = sim_targets(input_std, TARGET_TYPE_STANDARD, nscience) - nsky = sim_targets(input_sky, TARGET_TYPE_SKY, (nscience + nstd)) + input_suppsky = os.path.join(test_dir, "suppsky.fits") + tgoff = 0 + nscience = sim_targets( + input_mtl, + TARGET_TYPE_SCIENCE, + tgoff, + density=self.density_science + ) + tgoff += nscience + nstd = sim_targets( + input_std, + TARGET_TYPE_STANDARD, + tgoff, + density=self.density_standards + ) + tgoff += nstd + nsky = sim_targets( + input_sky, + TARGET_TYPE_SKY, + tgoff, + density=self.density_sky + ) + tgoff += nsky + nsuppsky = sim_targets( + input_suppsky, + TARGET_TYPE_SUPPSKY, + tgoff, + density=self.density_suppsky + ) tfile = os.path.join(test_dir, "footprint.fits") sim_tiles(tfile) opts = { - "targets": [input_mtl, input_std, input_sky], + "targets": [input_mtl, input_std, input_sky, input_suppsky], "dir": test_dir, "footprint": tfile, "standards_per_petal": 10, @@ -360,9 +456,11 @@ def test_cli(self): args = parse_assign(optlist) run_assign_full(args) + # plotpetals = "0" + plotpetals = "0,1,2,3,4,5,6,7,8,9" opts = { "dir": test_dir, - "petals": "0", + "petals": plotpetals, "serial": True } optlist = option_list(opts) @@ -389,7 +487,9 @@ def test_cli(self): for tile, props in qadata.items(): self.assertEqual(4500, props["assign_science"]) self.assertEqual(100, props["assign_std"]) - self.assertEqual(400, props["assign_sky"]) + self.assertEqual(400, ( + props["assign_sky"] + props["assign_suppsky"] + )) return def test_fieldrot(self): @@ -398,9 +498,35 @@ def test_fieldrot(self): input_mtl = os.path.join(test_dir, "mtl.fits") input_std = os.path.join(test_dir, "standards.fits") input_sky = os.path.join(test_dir, "sky.fits") - nscience = sim_targets(input_mtl, TARGET_TYPE_SCIENCE, 0) - nstd = sim_targets(input_std, TARGET_TYPE_STANDARD, nscience) - nsky = sim_targets(input_sky, TARGET_TYPE_SKY, (nscience + nstd)) + input_suppsky = os.path.join(test_dir, "suppsky.fits") + tgoff = 0 + nscience = sim_targets( + input_mtl, + TARGET_TYPE_SCIENCE, + tgoff, + density=self.density_science + ) + tgoff += nscience + nstd = sim_targets( + input_std, + TARGET_TYPE_STANDARD, + tgoff, + density=self.density_standards + ) + tgoff += nstd + nsky = sim_targets( + input_sky, + TARGET_TYPE_SKY, + tgoff, + density=self.density_sky + ) + tgoff += nsky + nsuppsky = sim_targets( + input_suppsky, + TARGET_TYPE_SUPPSKY, + tgoff, + density=self.density_suppsky + ) # Simulate the tiles tfile = os.path.join(test_dir, "footprint.fits") @@ -420,6 +546,7 @@ def test_fieldrot(self): load_target_file(tgs, input_mtl) load_target_file(tgs, input_std) load_target_file(tgs, input_sky) + load_target_file(tgs, input_suppsky) # Create a hierarchical triangle mesh lookup of the targets # positions diff --git a/py/fiberassign/test/test_targets.py b/py/fiberassign/test/test_targets.py index 7ed438b5..d5c6c81e 100644 --- a/py/fiberassign/test/test_targets.py +++ b/py/fiberassign/test/test_targets.py @@ -15,6 +15,7 @@ from fiberassign.tiles import load_tiles from fiberassign.targets import (TARGET_TYPE_SCIENCE, TARGET_TYPE_SKY, + TARGET_TYPE_SUPPSKY, TARGET_TYPE_STANDARD, load_target_file, desi_target_type, default_main_sciencemask, default_main_skymask, default_main_stdmask, @@ -40,14 +41,21 @@ def test_available(self): input_mtl = os.path.join(test_dir, "mtl.fits") input_std = os.path.join(test_dir, "standards.fits") input_sky = os.path.join(test_dir, "sky.fits") - nscience = sim_targets(input_mtl, TARGET_TYPE_SCIENCE, 0) - nstd = sim_targets(input_std, TARGET_TYPE_STANDARD, nscience) - nsky = sim_targets(input_sky, TARGET_TYPE_SKY, (nscience + nstd)) + input_suppsky = os.path.join(test_dir, "suppsky.fits") + tgoff = 0 + nscience = sim_targets(input_mtl, TARGET_TYPE_SCIENCE, tgoff) + tgoff += nscience + nstd = sim_targets(input_std, TARGET_TYPE_STANDARD, tgoff) + tgoff += nstd + nsky = sim_targets(input_sky, TARGET_TYPE_SKY, tgoff) + tgoff += nsky + nsuppsky = sim_targets(input_suppsky, TARGET_TYPE_SUPPSKY, tgoff) tgs = Targets() load_target_file(tgs, input_mtl) load_target_file(tgs, input_std) load_target_file(tgs, input_sky) + load_target_file(tgs, input_suppsky) print(tgs) # Create a hierarchical triangle mesh lookup of the targets positions @@ -77,12 +85,14 @@ def test_target_type(self): desi_mask["ELG"].mask, desi_mask["STD_FAINT"].mask, desi_mask["SKY"].mask, + desi_mask["SUPP_SKY"].mask, desi_mask["IN_BRIGHT_OBJECT"].mask, ] fbatype = np.array([ TARGET_TYPE_SCIENCE, TARGET_TYPE_STANDARD, TARGET_TYPE_SKY, + TARGET_TYPE_SUPPSKY, 0 ]) result = desi_target_type( diff --git a/src/assign.cpp b/src/assign.cpp index b6149900..63b0841c 100644 --- a/src/assign.cpp +++ b/src/assign.cpp @@ -52,6 +52,7 @@ fba::Assignment::Assignment(fba::Targets::pshr tgs, tgtypes.push_back(TARGET_TYPE_SCIENCE); tgtypes.push_back(TARGET_TYPE_STANDARD); tgtypes.push_back(TARGET_TYPE_SKY); + tgtypes.push_back(TARGET_TYPE_SUPPSKY); tgtypes.push_back(TARGET_TYPE_SAFE); tile_target_xy.clear(); @@ -351,14 +352,23 @@ void fba::Assignment::assign_unused(uint8_t tgtype, int32_t max_per_petal, // The petal int32_t p = hw_->loc_petal[lid]; - if (nassign_petal[tgtype][tile_id][p] >= max_per_petal) { + // Special handling of supplemental sky targets. When + // enforcing per-petal limits on assignment counts of SUPP_SKY, + // we want to include the current assignment counts of regular + // SKY targets as well. + int32_t cur_petal = nassign_petal[tgtype][tile_id][p]; + if (tgtype == TARGET_TYPE_SUPPSKY) { + cur_petal += nassign_petal[TARGET_TYPE_SKY][tile_id][p]; + } + + if (cur_petal >= max_per_petal) { // Already have enough objects on this petal if (extra_log) { logmsg.str(""); logmsg << "assign unused " << tgstr << ": tile " << tile_id << ", petal " << p << " location " << lid << " (fiber " << fid << ")" << " has " - << nassign_petal[tgtype][tile_id][p] + << cur_petal << " (>= " << max_per_petal << ")"; logger.debug_tfg(tile_id, lid, -1, logmsg.str().c_str()); } @@ -841,15 +851,23 @@ void fba::Assignment::assign_force(uint8_t tgtype, int32_t required_per_petal, auto const & target_xy = tile_target_xy.at(tile_id); for (int32_t p = 0; p < npetal; ++p) { + // Special handling of supplemental sky targets. When + // enforcing per-petal limits on assignment counts of SUPP_SKY, + // we want to include the current assignment counts of regular + // SKY targets as well. + int32_t cur_petal = nassign_petal[tgtype][tile_id][p]; + if (tgtype == TARGET_TYPE_SUPPSKY) { + cur_petal += nassign_petal[TARGET_TYPE_SKY][tile_id][p]; + } if (extra_log) { logmsg.str(""); logmsg << "assign force " << tgstr << ": tile " << tile_id << ", petal " << p << " has " - << nassign_petal[tgtype][tile_id][p] + << cur_petal << " (require " << required_per_petal << ")"; logger.debug_tfg(tile_id, -1, -1, logmsg.str().c_str()); } - if (nassign_petal[tgtype][tile_id][p] >= required_per_petal) { + if (cur_petal >= required_per_petal) { // Already have enough objects on this petal continue; } @@ -1147,34 +1165,43 @@ void fba::Assignment::assign_force(uint8_t tgtype, int32_t required_per_petal, } } + cur_petal = nassign_petal[tgtype][tile_id][p]; + if (tgtype == TARGET_TYPE_SUPPSKY) { + cur_petal += nassign_petal[TARGET_TYPE_SKY][tile_id][p]; + } if (extra_log) { logmsg.str(""); logmsg << "assign force " << tgstr << ": tile " << tile_id << ", petal " << p << " now has " - << nassign_petal[tgtype][tile_id][p] + <= required_per_petal) { + if (cur_petal >= required_per_petal) { // Have enough objects on this petal break; } } - if (nassign_petal[tgtype][tile_id][p] - >= required_per_petal) { + cur_petal = nassign_petal[tgtype][tile_id][p]; + if (tgtype == TARGET_TYPE_SUPPSKY) { + cur_petal += nassign_petal[TARGET_TYPE_SKY][tile_id][p]; + } + if (cur_petal >= required_per_petal) { // Have enough objects on this petal break; } } - if (nassign_petal[tgtype][tile_id][p] - < required_per_petal) { + cur_petal = nassign_petal[tgtype][tile_id][p]; + if (tgtype == TARGET_TYPE_SUPPSKY) { + cur_petal += nassign_petal[TARGET_TYPE_SKY][tile_id][p]; + } + if (cur_petal < required_per_petal) { logmsg.str(""); logmsg << "assign force " << tgstr << ": tile " << tile_id << ", petal " << p << " could only assign " - << nassign_petal[tgtype][tile_id][p] + << cur_petal << " (require " << required_per_petal << "). Insufficient number of objects or too many collisions"; logger.warning(logmsg.str().c_str()); @@ -1676,7 +1703,7 @@ void fba::Assignment::assign_tileloc(fba::Hardware const * hw, // to update the counts for valid types of this object. static const std::vector target_types = { TARGET_TYPE_SCIENCE, TARGET_TYPE_STANDARD, - TARGET_TYPE_SKY, TARGET_TYPE_SAFE}; + TARGET_TYPE_SKY, TARGET_TYPE_SUPPSKY, TARGET_TYPE_SAFE}; for (auto const & tt : target_types) { if (tgobj.is_type(tt)) { nassign_tile.at(tt).at(tile)++; @@ -1750,7 +1777,7 @@ void fba::Assignment::unassign_tileloc(fba::Hardware const * hw, // to update the counts for valid types of this object. static const std::vector target_types = { TARGET_TYPE_SCIENCE, TARGET_TYPE_STANDARD, - TARGET_TYPE_SKY, TARGET_TYPE_SAFE}; + TARGET_TYPE_SKY, TARGET_TYPE_SUPPSKY, TARGET_TYPE_SAFE}; for (auto const & tt : target_types) { if (tgobj.is_type(tt)) { nassign_tile.at(tt).at(tile)--;