diff --git a/gwemopt/args.py b/gwemopt/args.py index 0d617d0..047f9b4 100644 --- a/gwemopt/args.py +++ b/gwemopt/args.py @@ -38,7 +38,6 @@ def parse_args(args): parser.add_argument("--doIterativeTiling", action="store_true", default=False) parser.add_argument("--doMinimalTiling", action="store_true", default=False) parser.add_argument("--doOverlappingScheduling", action="store_true", default=False) - parser.add_argument("--doPerturbativeTiling", action="store_true", default=False) parser.add_argument("--doOrderByObservability", action="store_true", default=False) parser.add_argument("--catalogDir", help="catalog directory", default=CATALOG_DIR) diff --git a/gwemopt/coverage.py b/gwemopt/coverage.py index 1245aad..9ef925d 100644 --- a/gwemopt/coverage.py +++ b/gwemopt/coverage.py @@ -19,7 +19,6 @@ eject_tiles, optimize_max_tiles, order_by_observability, - perturb_tiles, slice_galaxy_tiles, slice_map_tiles, slice_number_tiles, @@ -159,15 +158,6 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): params, tile_struct, combine_coverage_structs(coverage_structs) ) - if ( - params["doPerturbativeTiling"] - and (jj > 0) - and (not params["tilesType"] == "galaxy") - ): - tile_struct = perturb_tiles( - params, config_struct, telescope, map_struct_hold, tile_struct - ) - if params["doOverlappingScheduling"]: tile_struct = check_overlapping_tiles( params, tile_struct, combine_coverage_structs(coverage_structs) diff --git a/gwemopt/moc.py b/gwemopt/moc.py index 852331b..655dba2 100644 --- a/gwemopt/moc.py +++ b/gwemopt/moc.py @@ -17,6 +17,57 @@ from gwemopt.utils.pixels import get_region_moc +def construct_moc(params, config_struct, tesselation): + + if params["doParallel"]: + n_threads = params["Ncores"] + else: + n_threads = None + + if params["doChipGaps"]: + if telescope == "ZTF": + mocs = get_ztf_quadrant_moc( + tesselation[:, 1], tesselation[:, 2], n_threads=n_threads + ) + elif telescope == "DECam": + mocs = get_decam_quadrant_moc( + tesselation[:, 1], tesselation[:, 2], n_threads=n_threads + ) + else: + raise ValueError("Chip gaps only available for DECam and ZTF") + + else: + if config_struct["FOV_type"] == "circle": + mocs = MOC.from_cones( + lon=tesselation[:, 1] * u.deg, + lat=tesselation[:, 2] * u.deg, + radius=config_struct["FOV"] * u.deg, + max_depth=np.uint8(10), + n_threads=n_threads, + ) + elif config_struct["FOV_type"] == "square": + mocs = MOC.from_boxes( + lon=tesselation[:, 1] * u.deg, + lat=tesselation[:, 2] * u.deg, + a=config_struct["FOV"] * u.deg, + b=config_struct["FOV"] * u.deg, + angle=0 * u.deg, + max_depth=np.uint8(10), + n_threads=n_threads, + ) + elif config_struct["FOV_type"] == "region": + region_file = os.path.join(CONFIG_DIR, config_struct["FOV"]) + region = regions.Regions.read(region_file, format="ds9") + mocs = get_region_moc( + tesselation[:, 1], tesselation[:, 2], region, n_threads=n_threads + ) + + return { + tess[0].astype(int): {"ra": tess[1], "dec": tess[2], "moc": mocs[ii]} + for ii, tess in enumerate(tesselation) + } + + def create_moc(params, map_struct=None): nside = params["nside"] @@ -24,7 +75,6 @@ def create_moc(params, map_struct=None): for telescope in params["telescopes"]: config_struct = params["config"][telescope] tesselation = config_struct["tesselation"] - moc_struct = {} if (telescope == "ZTF") and params["doUsePrimary"]: idx = np.where(tesselation[:, 0] <= 880)[0] @@ -33,58 +83,7 @@ def create_moc(params, map_struct=None): idx = np.where(tesselation[:, 0] >= 1000)[0] tesselation = tesselation[idx, :] - if params["doChipGaps"]: - if params["doParallel"]: - n_threads = params["Ncores"] - else: - n_threads = None - if telescope == "ZTF": - mocs = get_ztf_quadrant_moc( - tesselation[:, 1], tesselation[:, 2], n_threads=n_threads - ) - elif telescope == "DECam": - mocs = get_decam_quadrant_moc( - tesselation[:, 1], tesselation[:, 2], n_threads=n_threads - ) - else: - raise ValueError("Chip gaps only available for DECam and ZTF") - - else: - if params["doParallel"]: - n_threads = params["Ncores"] - else: - n_threads = None - - if config_struct["FOV_type"] == "circle": - mocs = MOC.from_cones( - lon=tesselation[:, 1] * u.deg, - lat=tesselation[:, 2] * u.deg, - radius=config_struct["FOV"] * u.deg, - max_depth=np.uint8(10), - n_threads=n_threads, - ) - elif config_struct["FOV_type"] == "square": - mocs = MOC.from_boxes( - lon=tesselation[:, 1] * u.deg, - lat=tesselation[:, 2] * u.deg, - a=config_struct["FOV"] * u.deg, - b=config_struct["FOV"] * u.deg, - angle=0 * u.deg, - max_depth=np.uint8(10), - n_threads=n_threads, - ) - elif config_struct["FOV_type"] == "region": - region_file = os.path.join(CONFIG_DIR, config_struct["FOV"]) - region = regions.Regions.read(region_file, format="ds9") - mocs = get_region_moc( - tesselation[:, 1], tesselation[:, 2], region, n_threads=n_threads - ) - - moc_struct = { - tess[0].astype(int): {"ra": tess[1], "dec": tess[2], "moc": mocs[ii]} - for ii, tess in enumerate(tesselation) - } - + moc_struct = construct_moc(params, config_struct, tesselation) if map_struct is not None: moc_keep = map_struct["moc_keep"] else: diff --git a/gwemopt/tiles.py b/gwemopt/tiles.py index cdfe9e6..8301dcb 100644 --- a/gwemopt/tiles.py +++ b/gwemopt/tiles.py @@ -491,88 +491,6 @@ def order_by_observability(params, tile_structs): params["telescopes"] = [params["telescopes"][ii] for ii in idx] -def perturb_tiles(params, config_struct, telescope, map_struct, tile_struct): - map_struct_hold = copy.deepcopy(map_struct) - ipix_keep = map_struct_hold["ipix_keep"] - nside = params["nside"] - - if config_struct["FOV_type"] == "square": - width = config_struct["FOV"] * 0.5 - elif config_struct["FOV_type"] == "circle": - width = config_struct["FOV"] * 1.0 - - moc_struct = {} - keys = list(tile_struct.keys()) - for ii, key in enumerate(keys): - if tile_struct[key]["prob"] == 0.0: - continue - - if np.mod(ii, 100) == 0: - print("Optimizing tile %d/%d" % (ii, len(keys))) - - x0 = [tile_struct[key]["ra"], tile_struct[key]["dec"]] - FOV = config_struct["FOV"] - bounds = [ - [tile_struct[key]["ra"] - width, tile_struct[key]["ra"] + width], - [tile_struct[key]["dec"] - width, tile_struct[key]["dec"] + width], - ] - - ras = np.linspace( - tile_struct[key]["ra"] - width, tile_struct[key]["ra"] + width, 5 - ) - decs = np.linspace( - tile_struct[key]["dec"] - width, tile_struct[key]["dec"] + width, 5 - ) - RAs, DECs = np.meshgrid(ras, decs) - ras, decs = RAs.flatten(), DECs.flatten() - - vals = [] - for ra, dec in zip(ras, decs): - if np.abs(dec) > 90: - vals.append(0) - continue - - moc_struct_temp = gwemopt.moc.Fov2Moc( - params, config_struct, telescope, ra, dec, nside - ) - idx = np.where(map_struct_hold["prob"][moc_struct_temp["ipix"]] == -1)[0] - idx = np.setdiff1d(idx, ipix_keep) - if len(map_struct_hold["prob"][moc_struct_temp["ipix"]]) == 0: - rat = 0.0 - else: - rat = float(len(idx)) / float( - len(map_struct_hold["prob"][moc_struct_temp["ipix"]]) - ) - if rat > params["maximumOverlap"]: - val = 0.0 - else: - ipix = moc_struct_temp["ipix"] - if len(ipix) == 0: - val = 0.0 - else: - vals_to_sum = map_struct_hold["prob"][ipix] - vals_to_sum[vals_to_sum < 0] = 0 - val = np.sum(vals_to_sum) - vals.append(val) - idx = np.argmax(vals) - ra, dec = ras[idx], decs[idx] - moc_struct[key] = gwemopt.moc.Fov2Moc( - params, config_struct, telescope, ra, dec, nside - ) - - map_struct_hold["prob"][moc_struct[key]["ipix"]] = -1 - ipix_keep = np.setdiff1d(ipix_keep, moc_struct[key]["ipix"]) - - tile_struct = powerlaw_tiles_struct( - params, config_struct, telescope, map_struct, moc_struct - ) - tile_struct = gwemopt.segments.get_segments_tiles( - params, config_struct, tile_struct - ) - - return tile_struct - - def schedule_alternating( params, config_struct, @@ -833,12 +751,16 @@ def galaxy(params, map_struct, catalog_struct: pd.DataFrame): # redefine catalog_struct catalog_struct_new = pd.DataFrame(new_cat) - moc_struct = {} + tesselation = np.vstack( + ( + np.arange(len(catalog_struct_new["ra"])), + catalog_struct_new["ra"], + catalog_struct_new["dec"], + ) + ).T + moc_struct = gwemopt.moc.construct_moc(params, config_struct, tesselation) cnt = 0 for _, row in catalog_struct_new.iterrows(): - moc_struct[cnt] = gwemopt.moc.Fov2Moc( - params, config_struct, telescope, row["ra"], row["dec"], nside - ) moc_struct[cnt]["galaxies"] = row["galaxies"] cnt = cnt + 1