Skip to content

Commit

Permalink
use new functionality for galaxies as well
Browse files Browse the repository at this point in the history
  • Loading branch information
@mcoughlin
mcoughlin committed Jun 25, 2024
1 parent 288baa5 commit e1b3845
Show file tree
Hide file tree
Showing 4 changed files with 60 additions and 150 deletions.
1 change: 0 additions & 1 deletion gwemopt/args.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
10 changes: 0 additions & 10 deletions gwemopt/coverage.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@
eject_tiles,
optimize_max_tiles,
order_by_observability,
perturb_tiles,
slice_galaxy_tiles,
slice_map_tiles,
slice_number_tiles,
Expand Down Expand Up @@ -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)
Expand Down
105 changes: 52 additions & 53 deletions gwemopt/moc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,14 +17,64 @@
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"]

moc_structs = {}
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]
Expand All @@ -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:
Expand Down
94 changes: 8 additions & 86 deletions gwemopt/tiles.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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

Expand Down

0 comments on commit e1b3845

Please sign in to comment.