From bbd534b3fd509c5e5269c09d3183138af411928e Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Tue, 28 May 2024 09:54:58 +0200 Subject: [PATCH] wip --- gwemopt/args.py | 7 +- gwemopt/coverage.py | 44 +----- gwemopt/io/__init__.py | 5 +- gwemopt/io/export_tiles.py | 10 -- gwemopt/io/schedule.py | 47 +++--- gwemopt/io/skymap.py | 223 ++++++++++++--------------- gwemopt/moc.py | 18 +-- gwemopt/params.py | 5 +- gwemopt/plotting/plot_coverage.py | 159 +++----------------- gwemopt/plotting/plot_skymap.py | 64 ++------ gwemopt/plotting/plot_tiles.py | 76 +++------- gwemopt/plotting/style.py | 28 ++-- gwemopt/scheduler.py | 172 +-------------------- gwemopt/tiles.py | 241 +++++++----------------------- gwemopt/utils/pixels.py | 52 ++----- 15 files changed, 284 insertions(+), 867 deletions(-) diff --git a/gwemopt/args.py b/gwemopt/args.py index 6aceb9b1..ee5af211 100644 --- a/gwemopt/args.py +++ b/gwemopt/args.py @@ -83,7 +83,6 @@ def parse_args(args): parser.add_argument("--Ndet", default=1, type=int) parser.add_argument("--nside", default=256, type=int) - parser.add_argument("--DScale", default=1.0, type=float) parser.add_argument("--Tobs", default="0.0,1.0") parser.add_argument("--mindiff", default=0.0, type=float) @@ -111,8 +110,6 @@ def parse_args(args): parser.add_argument("--doSingleExposure", action="store_true", default=False) parser.add_argument("--filters", default="r,g,r") parser.add_argument("--doAlternatingFilters", action="store_true", default=False) - parser.add_argument("--doRASlices", action="store_true", default=False) - parser.add_argument("--nside_down", default=2, type=int) parser.add_argument("--max_filter_sets", default=4, type=int) parser.add_argument("--iterativeOverlap", default=0.0, type=float) parser.add_argument("--maximumOverlap", default=1.0, type=float) @@ -137,10 +134,8 @@ def parse_args(args): parser.add_argument("--doBlocks", action="store_true", default=False) parser.add_argument("--Nblocks", default=4, type=int) - parser.add_argument("--doRASlice", action="store_true", default=False) - parser.add_argument("--raslice", default="0.0,24.0") - parser.add_argument("--absmag", default=-15.0, type=float) + parser.add_argument("--raslice", default=None) parser.add_argument( "--galactic_limit", help="Galactic limit.", default=0.0, type=float diff --git a/gwemopt/coverage.py b/gwemopt/coverage.py index d3fa0b59..985c88ec 100644 --- a/gwemopt/coverage.py +++ b/gwemopt/coverage.py @@ -30,10 +30,7 @@ def combine_coverage_structs(coverage_structs): coverage_struct_combined = {} coverage_struct_combined["data"] = np.empty((0, 8)) coverage_struct_combined["filters"] = np.empty((0, 1)) - coverage_struct_combined["ipix"] = [] - coverage_struct_combined["patch"] = [] - coverage_struct_combined["FOV"] = np.empty((0, 1)) - coverage_struct_combined["area"] = np.empty((0, 1)) + coverage_struct_combined["moc"] = [] coverage_struct_combined["telescope"] = np.empty((0, 1)) coverage_struct_combined["galaxies"] = [] coverage_struct_combined["exposureused"] = [] @@ -45,17 +42,8 @@ def combine_coverage_structs(coverage_structs): coverage_struct_combined["filters"] = np.append( coverage_struct_combined["filters"], coverage_struct["filters"] ) - coverage_struct_combined["ipix"] = ( - coverage_struct_combined["ipix"] + coverage_struct["ipix"] - ) - coverage_struct_combined["patch"] = ( - coverage_struct_combined["patch"] + coverage_struct["patch"] - ) - coverage_struct_combined["FOV"] = np.append( - coverage_struct_combined["FOV"], coverage_struct["FOV"] - ) - coverage_struct_combined["area"] = np.append( - coverage_struct_combined["area"], coverage_struct["area"] + coverage_struct_combined["moc"] = ( + coverage_struct_combined["moc"] + coverage_struct["moc"] ) coverage_struct_combined["telescope"] = np.append( coverage_struct_combined["telescope"], coverage_struct["telescope"] @@ -173,7 +161,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): coverage_structs = [] n_scope = 0 - full_prob_map = map_struct["prob"] + full_prob_map = map_struct["skymap"] for jj, telescope in enumerate(params["telescopes"]): if params["splitType"] is not None: @@ -233,16 +221,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): previous_coverage_struct, ) - if params["doRASlices"]: - coverage_struct = gwemopt.scheduler.schedule_ra_splits( - params, - config_struct, - map_struct_hold, - tile_struct, - telescope, - previous_coverage_struct, - ) - elif params["doBalanceExposure"]: + if params["doBalanceExposure"]: if ( params["max_nb_tiles"] is None ): # optimize max tiles (iff max tiles not already specified) @@ -435,7 +414,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): if params["doIterativeTiling"]: map_struct_hold = slice_map_tiles(params, map_struct_hold, coverage_struct) - map_struct["prob"] = full_prob_map + map_struct["skymap"] = full_prob_map return tile_structs, combine_coverage_structs(coverage_structs) @@ -508,16 +487,7 @@ def absmag(params, map_struct, tile_structs, previous_coverage_struct=None): tile_struct, previous_coverage_struct, ) - if params["doRASlices"]: - coverage_struct = gwemopt.scheduler.schedule_ra_splits( - params, - config_struct, - map_struct_hold, - tile_struct, - telescope, - previous_coverage_struct, - ) - elif params["doBalanceExposure"]: + if params["doBalanceExposure"]: if params["max_nb_tiles"] is None: # optimize max tiles (iff max tiles not already specified) optimized_max, coverage_struct, tile_struct = optimize_max_tiles( diff --git a/gwemopt/io/__init__.py b/gwemopt/io/__init__.py index 6a59ac8e..758e542d 100644 --- a/gwemopt/io/__init__.py +++ b/gwemopt/io/__init__.py @@ -1,6 +1,3 @@ -from gwemopt.io.export_tiles import ( - export_tiles_coverage_hist, - export_tiles_coverage_int, -) +from gwemopt.io.export_tiles import export_tiles_coverage_int from gwemopt.io.schedule import export_schedule_xml, read_schedule, summary from gwemopt.io.skymap import get_skymap, read_skymap diff --git a/gwemopt/io/export_tiles.py b/gwemopt/io/export_tiles.py index e6926c03..50112696 100644 --- a/gwemopt/io/export_tiles.py +++ b/gwemopt/io/export_tiles.py @@ -1,13 +1,3 @@ -def export_tiles_coverage_hist(params, diffs): - """ - Export the histogram of the tiles coverage. - """ - filename = params["outputDir"].joinpath("tiles_coverage_hist.dat") - with open(filename, "w") as fid: - for ii in range(len(diffs)): - fid.write("%.10f\n" % diffs[ii]) - - def export_tiles_coverage_int(filename, tts, cum_probs, cum_areas): with open(filename, "w") as fid: for ii in range(len(tts)): diff --git a/gwemopt/io/schedule.py b/gwemopt/io/schedule.py index a6745d75..1345681e 100644 --- a/gwemopt/io/schedule.py +++ b/gwemopt/io/schedule.py @@ -151,13 +151,10 @@ def export_schedule_xml(xmlfile, map_struct, coverage_struct, config_struct): Field(name="priority", ucd="", unit="", dataType="int", Description=[""]) ) table_field = utilityTable(table) - table_field.blankTable(len(coverage_struct["ipix"])) + table_field.blankTable(len(coverage_struct["moc"])) - for ii in range(len(coverage_struct["ipix"])): + for ii in range(len(coverage_struct["moc"])): data = coverage_struct["data"][ii, :] - ipix = coverage_struct["ipix"][ii] - - prob = np.sum(map_struct["prob"][ipix]) ra, dec = data[0], data[1] observ_time, exposure_time, field_id, prob, airmass = ( @@ -198,15 +195,6 @@ def export_schedule_xml(xmlfile, map_struct, coverage_struct, config_struct): def summary(params, map_struct, coverage_struct, catalog_struct=None): - idx50 = len(map_struct["cumprob"]) - np.argmin(np.abs(map_struct["cumprob"] - 0.50)) - idx90 = len(map_struct["cumprob"]) - np.argmin(np.abs(map_struct["cumprob"] - 0.90)) - - mapfile = params["outputDir"].joinpath("map.dat") - with open(mapfile, "w") as fid: - fid.write( - "%.5f %.5f\n" - % (map_struct["pixarea_deg2"] * idx50, map_struct["pixarea_deg2"] * idx90) - ) filts = list(set(coverage_struct["filters"])) for jj, telescope in enumerate(params["telescopes"]): @@ -225,15 +213,12 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): totexp = 0 with open(schedulefile, "w") as fid: - for ii in range(len(coverage_struct["ipix"])): + for ii in range(len(coverage_struct["filters"])): if not telescope == coverage_struct["telescope"][ii]: continue data = coverage_struct["data"][ii, :] filt = coverage_struct["filters"][ii] - ipix = coverage_struct["ipix"][ii] - - prob = np.sum(map_struct["prob"][ipix]) ra, dec = data[0], data[1] observ_time, mag, exposure_time, field_id, prob, airmass = ( @@ -307,6 +292,8 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): fid.write("\n") summaryfile = params["outputDir"].joinpath("summary.dat") + cummoc = None + with open(summaryfile, "w") as fid: gpstime = params["gpstime"] event_mjd = Time(gpstime, format="gps", scale="utc").mjd @@ -322,26 +309,36 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): if params["tilesType"] == "galaxy": galaxies = np.empty((0, 2)) - for ii in range(len(coverage_struct["ipix"])): + for ii in range(len(coverage_struct["filters"])): data = coverage_struct["data"][ii, :] filt = coverage_struct["filters"][ii] - ipix = coverage_struct["ipix"][ii] + moc = coverage_struct["moc"][ii] - prob = np.sum(map_struct["prob"][ipix]) + ra, dec = data[0], data[1] + observ_time, mag, exposure_time, field_id, prob, airmass = ( + data[2], + data[3], + data[4], + data[5], + data[6], + data[7], + ) if data[2] > event_mjd + tt: continue - ipixs = np.append(ipixs, ipix) - ipixs = np.unique(ipixs).astype(int) - cum_prob = np.sum(map_struct["prob"][ipixs]) + if cummoc is None: + cummoc = moc + else: + cummoc = cummoc + moc + cum_prob = 0 ## FIXME if params["tilesType"] == "galaxy": galaxies = np.append(galaxies, coverage_struct["galaxies"][ii]) galaxies = np.unique(galaxies).astype(int) cum_prob = np.sum(catalog_struct[params["galaxy_grade"]][galaxies]) - cum_area = len(ipixs) * map_struct["pixarea_deg2"] + cum_area = 0 ## FIXME mjds.append(data[2]) mjds_floor.append(int(np.floor(data[2]))) diff --git a/gwemopt/io/skymap.py b/gwemopt/io/skymap.py index 9bf59e98..c95c6af5 100644 --- a/gwemopt/io/skymap.py +++ b/gwemopt/io/skymap.py @@ -5,8 +5,10 @@ import os from pathlib import Path +import astropy_healpix as ah import healpy as hp import ligo.skymap.distance as ligodist +import ligo.skymap.plot import lxml.etree import numpy as np import requests @@ -15,9 +17,11 @@ from astropy.io import fits from astropy.time import Time from ligo.gracedb.rest import GraceDb +from ligo.skymap import moc from ligo.skymap.bayestar import rasterize from ligo.skymap.io import read_sky_map -from ligo.skymap.moc import uniq2nest +from matplotlib import pyplot as plt +from mocpy import MOC from scipy.interpolate import PchipInterpolator from scipy.stats import norm @@ -137,14 +141,14 @@ def read_inclination(skymap, params, map_struct): theta = 0.5 * np.pi - np.deg2rad(dec) # make use of the maP nside maP_idx = np.argmax(skymap["PROBDENSITY"]) - order, _ = uniq2nest(skymap[maP_idx]["UNIQ"]) + order, _ = moc.uniq2nest(skymap[maP_idx]["UNIQ"]) nside = hp.order2nside(order) # get the nested idx for the given sky location nest_idx = hp.ang2pix(nside, theta, phi, nest=True) # find the row with the closest nested index nest_idxs = [] for row in skymap: - order_per_row, nest_idx_per_row = uniq2nest(row["UNIQ"]) + order_per_row, nest_idx_per_row = moc.uniq2nest(row["UNIQ"]) if order_per_row == order: nest_idxs.append(nest_idx_per_row) else: @@ -156,7 +160,7 @@ def read_inclination(skymap, params, map_struct): maP_idx = np.argmax(skymap["PROBDENSITY"]) uniq_idx = skymap[maP_idx]["UNIQ"] # convert to nested indexing and find the location of that index - order, nest_idx = uniq2nest(uniq_idx) + order, nest_idx = moc.uniq2nest(uniq_idx) nside = hp.order2nside(order) theta, phi = hp.pix2ang(nside, int(nest_idx), nest=True) # convert theta and phi to ra and dec @@ -258,151 +262,114 @@ def read_skymap(params, map_struct=None): if "do_3d" not in params: params["do_3d"] = is_3d - header = [] - if map_struct is None: - map_struct = {} - - filename = params["skymap"] - - if params["do_3d"]: - try: - healpix_data, header = hp.read_map( - filename, field=(0, 1, 2, 3), h=True - ) - except: - skymap = read_sky_map(filename, moc=True, distances=True) - order = hp.nside2order(params["nside"]) - - # for colname in skymap.colnames: - # if colname.startswith('PROB'): - # newname = colname.replace('PROB', 'PROBDENSITY') - # skymap.rename_column(colname, newname) - # skymap[newname] *= len(skymap) / (4 * np.pi) - # skymap[newname].unit = u.steradian ** -1 - - if "PROBDENSITY_SAMPLES" in skymap.columns: - if params["inclination"]: - map_struct = read_inclination(skymap, params, map_struct) - - skymap.remove_columns( - [ - f"{name}_SAMPLES" - for name in [ - "PROBDENSITY", - "DISTMU", - "DISTSIGMA", - "DISTNORM", - ] - ] - ) - - t = rasterize(skymap) - result = t["PROB"], t["DISTMU"], t["DISTSIGMA"], t["DISTNORM"] - healpix_data = hp.reorder(result, "NESTED", "RING") - - distmu_data = healpix_data[1] - distsigma_data = healpix_data[2] - prob_data = healpix_data[0] - norm_data = healpix_data[3] - - map_struct["distmu"] = distmu_data / params["DScale"] - map_struct["distsigma"] = distsigma_data / params["DScale"] - map_struct["prob"] = prob_data - map_struct["distnorm"] = norm_data + map_struct = {} - else: - prob_data, header = hp.read_map( - filename, field=0, verbose=False, h=True - ) - prob_data = prob_data / np.sum(prob_data) - - map_struct["prob"] = prob_data - - for j in range(len(header)): - if header[j][0] == "DATE": - map_struct["trigtime"] = header[j][1] + filename = params["skymap"] - natural_nside = hp.pixelfunc.get_nside(map_struct["prob"]) - nside = params["nside"] + if params["do_3d"]: + skymap = read_sky_map(filename, moc=True, distances=True) - print("natural_nside =", natural_nside) - print("nside =", nside) + if "PROBDENSITY_SAMPLES" in skymap.columns: + if params["inclination"]: + map_struct = read_inclination(skymap, params, map_struct) - if not params["do_3d"]: - map_struct["prob"] = hp.ud_grade(map_struct["prob"], nside, power=-2) - - if params["do_3d"]: - if natural_nside != nside: - map_struct["prob"] = hp.pixelfunc.ud_grade( - map_struct["prob"], nside, power=-2 - ) - map_struct["distmu"] = hp.pixelfunc.ud_grade(map_struct["distmu"], nside) - map_struct["distsigma"] = hp.pixelfunc.ud_grade( - map_struct["distsigma"], nside - ) - map_struct["distnorm"] = hp.pixelfunc.ud_grade( - map_struct["distnorm"], nside - ) + skymap.remove_columns( + [ + f"{name}_SAMPLES" + for name in [ + "PROBDENSITY", + "DISTMU", + "DISTSIGMA", + "DISTNORM", + ] + ] + ) - map_struct["distmu"][map_struct["distmu"] < -1e30] = np.inf + map_struct["skymap"] = skymap + else: + skymap = read_sky_map(filename, moc=True, distances=False) + map_struct["skymap"] = skymap + + level, ipix = ah.uniq_to_level_ipix(map_struct["skymap"]["UNIQ"]) + LEVEL = MOC.MAX_ORDER + shift = 2 * (LEVEL - level) + hpx = np.array(np.vstack([ipix << shift, (ipix + 1) << shift]), dtype=np.uint64).T + map_struct["skymap"]["HPX"] = hpx + nside = ah.level_to_nside(level) + pixel_area = ah.nside_to_pixel_area(ah.level_to_nside(level)) + ra, dec = ah.healpix_to_lonlat(ipix, nside, order="nested") + map_struct["skymap"]["ra"] = ra.deg + map_struct["skymap"]["dec"] = dec.deg + + if params["raslice"] is not None: + ra_low, ra_high = params["raslice"][0], params["raslice"][1] + if ra_low <= ra_high: + ipix = np.where( + (ra_high * 360.0 / 24.0 < ra.deg) | (ra_low * 360.0 / 24.0 > ra.deg) + )[0] + else: + ipix = np.where( + (ra_high * 360.0 / 24.0 < ra.deg) & (ra_low * 360.0 / 24.0 > ra.deg) + )[0] + map_struct["skymap"]["PROBDENSITY"][ipix] = 0.0 - nside_down = 32 + if params["galactic_limit"] > 0.0: + coords = SkyCoord(ra=ra, dec=dec) + ipix = np.where(np.abs(coords.galactic.b.deg) <= params["galactic_limit"])[0] + map_struct["skymap"]["PROBDENSITY"][ipix] = 0.0 - distmu_down = hp.pixelfunc.ud_grade(map_struct["distmu"], nside_down) + map_struct["skymap_raster"] = rasterize(map_struct["skymap"]) + if "DISTMU" in map_struct["skymap_raster"].columns: ( - map_struct["distmed"], - map_struct["diststd"], + map_struct["skymap_raster"]["DISTMEAN"], + map_struct["skymap_raster"]["DISTSTD"], mom_norm, ) = ligodist.parameters_to_moments( - map_struct["distmu"], map_struct["distsigma"] + map_struct["skymap_raster"]["DISTMU"], + map_struct["skymap_raster"]["DISTSIGMA"], ) - distmu_down[distmu_down < -1e30] = np.inf + extra_header = [ + ("PIXTYPE", "HEALPIX", "HEALPIX pixelisation"), + ("ORDERING", "NESTED", "Pixel ordering scheme: RING, NESTED, or NUNIQ"), + ("COORDSYS", "C", "Ecliptic, Galactic or Celestial (equatorial)"), + ( + "MOCORDER", + moc.uniq2order(map_struct["skymap"]["UNIQ"].max()), + "MOC resolution (best order)", + ), + ("INDXSCHM", "EXPLICIT", "Indexing: IMPLICIT or EXPLICIT"), + ] - map_struct["distmed"] = hp.ud_grade(map_struct["distmed"], nside, power=-2) - map_struct["diststd"] = hp.ud_grade(map_struct["diststd"], nside, power=-2) + hdu = fits.table_to_hdu(map_struct["skymap_raster"]) + hdu.header.extend(extra_header) + map_struct["hdu"] = hdu - npix = hp.nside2npix(nside) - theta, phi = hp.pix2ang(nside, np.arange(npix)) - ra = np.rad2deg(phi) - dec = np.rad2deg(0.5 * np.pi - theta) + skymap_sorted = map_struct["skymap"].copy() + skymap_sorted.sort("PROBDENSITY") + skymap_sorted.reverse() - map_struct["ra"] = ra - map_struct["dec"] = dec + level, ipix = ah.uniq_to_level_ipix(skymap_sorted["UNIQ"]) + pixel_area = ah.nside_to_pixel_area(ah.level_to_nside(level)) - if params["doRASlice"]: - ra_low, ra_high = params["raslice"][0], params["raslice"][1] - if ra_low <= ra_high: - ipix = np.where( - (ra_high * 360.0 / 24.0 < ra) | (ra_low * 360.0 / 24.0 > ra) - )[0] - else: - ipix = np.where( - (ra_high * 360.0 / 24.0 < ra) & (ra_low * 360.0 / 24.0 > ra) - )[0] - map_struct["prob"][ipix] = 0.0 - map_struct["prob"] = map_struct["prob"] / np.sum(map_struct["prob"]) + prob = pixel_area * skymap_sorted["PROBDENSITY"] + cumprob = np.cumsum(prob) - if params["galactic_limit"] > 0.0: - coords = SkyCoord(ra=ra * u.deg, dec=dec * u.deg) - ipix = np.where(np.abs(coords.galactic.b.deg) <= params["galactic_limit"])[0] - map_struct["prob"][ipix] = 0.0 - map_struct["prob"] = map_struct["prob"] / np.sum(map_struct["prob"]) + i = cumprob.searchsorted(params["iterativeOverlap"]) + skymap_keep = skymap_sorted[:i] - sort_idx = np.argsort(map_struct["prob"])[::-1] - csm = np.empty(len(map_struct["prob"])) - csm[sort_idx] = np.cumsum(map_struct["prob"][sort_idx]) + if params["iterativeOverlap"] > 0: - map_struct["cumprob"] = csm - map_struct["ipix_keep"] = np.where(csm <= params["iterativeOverlap"])[0] + tab = rasterize(skymap_keep) + hdu = fits.table_to_hdu(tab) + hdu.header.extend(extra_header) - pixarea = hp.nside2pixarea(nside) - pixarea_deg2 = hp.nside2pixarea(nside, degrees=True) + fig = plt.figure(figsize=(4, 4), dpi=100) + ax = plt.axes([0.05, 0.05, 0.9, 0.9], projection="astro mollweide") + ax.imshow_hpx(hdu, cmap="cylon") + fig.savefig("keep.pdf") - map_struct["nside"] = nside - map_struct["npix"] = npix - map_struct["pixarea"] = pixarea - map_struct["pixarea_deg2"] = pixarea_deg2 + map_struct["moc_keep"] = skymap_keep return params, map_struct diff --git a/gwemopt/moc.py b/gwemopt/moc.py index 796435f1..0ed13615 100644 --- a/gwemopt/moc.py +++ b/gwemopt/moc.py @@ -77,18 +77,19 @@ def create_moc(params, map_struct=None): ) if map_struct is not None: - ipix_keep = map_struct["ipix_keep"] + moc_keep = map_struct["moc_keep"] else: - ipix_keep = [] + moc_keep = None + if params["doMinimalTiling"]: moc_struct_new = copy.copy(moc_struct) if params["tilesType"] == "galaxy": tile_probs = gwemopt.tiles.compute_tiles_map( - params, moc_struct_new, prob, func="center", ipix_keep=ipix_keep + params, moc_struct_new, prob, func="center", moc_keep=moc_keep ) else: tile_probs = gwemopt.tiles.compute_tiles_map( - params, moc_struct_new, prob, func="np.sum(x)", ipix_keep=ipix_keep + params, moc_struct_new, prob, func="np.sum(x)", moc_keep=moc_keep ) keys = moc_struct.keys() @@ -130,7 +131,7 @@ def Fov2Moc(params, config_struct, telescope, ra_pointing, dec_pointing, nside): rotation = None if config_struct["FOV_type"] == "square": - ipix, radecs, patch, area = getSquarePixels( + moc = getSquarePixels( ra_pointing, dec_pointing, config_struct["FOV"], nside, rotation=rotation ) elif config_struct["FOV_type"] == "circle": @@ -158,11 +159,6 @@ def Fov2Moc(params, config_struct, telescope, ra_pointing, dec_pointing, nside): moc_struct["ra"] = ra_pointing moc_struct["dec"] = dec_pointing - moc_struct["ipix"] = ipix - moc_struct["corners"] = radecs - moc_struct["patch"] = patch - moc_struct["area"] = area - - moc_struct["moc"] = [] + moc_struct["moc"] = moc return moc_struct diff --git a/gwemopt/params.py b/gwemopt/params.py index 07211a2a..d18a093d 100644 --- a/gwemopt/params.py +++ b/gwemopt/params.py @@ -98,7 +98,10 @@ def params_struct(opts): params["treasuremap_status"] = opts.treasuremap_status.split(",") - params["raslice"] = np.array(opts.raslice.split(","), dtype=float) + if opts.raslice is not None: + params["raslice"] = np.array(opts.raslice.split(","), dtype=float) + else: + params["raslice"] = None params["unbalanced_tiles"] = None params["filters"] = opts.filters.split(",") diff --git a/gwemopt/plotting/plot_coverage.py b/gwemopt/plotting/plot_coverage.py index b888a1de..3aabddbb 100644 --- a/gwemopt/plotting/plot_coverage.py +++ b/gwemopt/plotting/plot_coverage.py @@ -6,56 +6,12 @@ from astropy.time import Time from matplotlib.pyplot import cm -from gwemopt.io import export_tiles_coverage_hist, export_tiles_coverage_int +from gwemopt.io import export_tiles_coverage_int from gwemopt.plotting.movie import make_movie from gwemopt.plotting.style import CBAR_BOOL, UNIT, add_edges, add_sun_moon, cmap from gwemopt.utils.geometry import angular_distance -def plot_mollweide_coverage(params, map_struct, coverage_struct): - """ - Plot the coverage in Mollweide projection. - """ - plot_name = params["outputDir"].joinpath("mollview_coverage.pdf") - plt.figure() - hp.mollview(map_struct["prob"], title="", unit=UNIT, cbar=CBAR_BOOL, cmap=cmap) - hp.projplot( - coverage_struct["data"][:, 0], - coverage_struct["data"][:, 1], - "wx", - lonlat=True, - coord="G", - ) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() - - -def plot_coverage(params, coverage_struct): - """ - Plot the coverage - """ - plot_name = params["outputDir"].joinpath("coverage.pdf") - plt.figure(figsize=(10, 8)) - for ii in range(len(coverage_struct["ipix"])): - data = coverage_struct["data"][ii, :] - filt = coverage_struct["filters"][ii] - - if filt == "g": - color = "g" - elif filt == "r": - color = "r" - else: - color = "k" - - plt.scatter(data[2], data[5], s=20, color=color) - - plt.xlabel("Time [MJD]") - plt.ylabel("Tile Number") - plt.savefig(plot_name, dpi=200) - plt.close() - - def plot_tiles_coverage(params, map_struct, coverage_struct, plot_sun_moon=False): """ Plot the tiles coverage in Mollweide projection. @@ -88,20 +44,6 @@ def plot_tiles_coverage(params, map_struct, coverage_struct, plot_sun_moon=False plt.close() -def plot_tiles_coverage_hist(params, diffs): - """ - Plot the histogram of the tiles coverage. - """ - plot_name = params["outputDir"].joinpath("tiles_coverage_hist.pdf") - plt.figure(figsize=(12, 8)) - bins = np.linspace(0.0, 24.0, 25) - plt.hist(24.0 * np.array(diffs), bins=bins) - plt.xlabel("Difference Between Observations [hours]") - plt.ylabel("Number of Observations") - plt.savefig(plot_name, dpi=200) - plt.close("all") - - def plot_tiles_coverage_int( params, map_struct, catalog_struct, coverage_struct, plot_sun_moon=False ): @@ -110,6 +52,9 @@ def plot_tiles_coverage_int( colors = cm.rainbow(np.linspace(0, 1, len(params["telescopes"]))) + hdu = map_struct["hdu"] + columns = [col.name for col in hdu.columns] + plot_name = params["outputDir"].joinpath("tiles_coverage_int.pdf") fig = plt.figure(figsize=(12, 8)) @@ -119,10 +64,7 @@ def plot_tiles_coverage_int( ax3 = ax2.twinx() # mirror them plt.axes(ax1) - hp.mollview( - map_struct["prob"], title="", unit=UNIT, cbar=CBAR_BOOL, cmap=cmap, hold=True - ) - add_edges() + ax1.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon") ax = plt.gca() if params["tilesType"] == "galaxy": @@ -136,30 +78,19 @@ def plot_tiles_coverage_int( color=color, ) else: - for ii in range(len(coverage_struct["ipix"])): - patch = coverage_struct["patch"][ii] - - idx = params["telescopes"].index(coverage_struct["telescope"][ii]) - - if patch == []: - continue - - if not type(patch) == list: - patch = [patch] - - for p in patch: - patch_cpy = copy.copy(p) - patch_cpy.axes = None - patch_cpy.figure = None - patch_cpy.set_transform(ax.transData) - patch_cpy.set_facecolor(colors[idx]) - - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch_cpy) + for ii in range(len(coverage_struct["moc"])): + moc = coverage_struct["moc"][ii] + data = coverage_struct["data"] + print(data) + moc.fill( + ax=ax, wcs=ax.wcs, alpha=data[0], fill=True, color="black", linewidth=1 + ) + moc.border(ax=ax, wcs=ax.wcs, alpha=1, color="black") idxs = np.argsort(coverage_struct["data"][:, 2]) plt.axes(ax2) for telescope, color in zip(params["telescopes"], colors): - ipixs = np.empty((0, 2)) + moc = None cum_prob = 0.0 tts, cum_probs, cum_areas = [], [], [] @@ -171,7 +102,7 @@ def plot_tiles_coverage_int( print("%s: %d/%d" % (telescope, jj, len(idxs))) data = coverage_struct["data"][ii, :] - ipix = coverage_struct["ipix"][ii] + m = coverage_struct["moc"][ii] if params["tilesType"] == "galaxy": galaxies = coverage_struct["galaxies"][ii] @@ -196,8 +127,10 @@ def plot_tiles_coverage_int( cum_galaxies = np.unique(cum_galaxies).astype(int) cum_area = len(cum_galaxies) else: - ipixs = np.append(ipixs, ipix) - ipixs = np.unique(ipixs).astype(int) + if moc is None: + moc = m + else: + moc = moc.union(m) cum_prob = np.sum(map_struct["prob"][ipixs]) cum_area = len(ipixs) * map_struct["pixarea_deg2"] @@ -380,7 +313,6 @@ def plot_coverage_scaled(params, map_struct, coverage_struct, plot_sun_moon, max def make_coverage_plots( params, map_struct, coverage_struct, catalog_struct=None, plot_sun_moon: bool = True ): - plot_mollweide_coverage(params, map_struct, coverage_struct) idx = np.isfinite(coverage_struct["data"][:, 4]) if not idx.size: @@ -388,59 +320,6 @@ def make_coverage_plots( max_time = np.max(coverage_struct["data"][idx, 4]) - plot_coverage(params, coverage_struct) - - plot_tiles_coverage( - params, map_struct, coverage_struct, plot_sun_moon=plot_sun_moon - ) - - diffs = [] - if params["tilesType"] == "galaxy": - coverage_ras = coverage_struct["data"][:, 0] - coverage_decs = coverage_struct["data"][:, 1] - - for ii in range(len(coverage_ras) - 1): - current_ra, current_dec = coverage_ras[ii], coverage_decs[ii] - - dist = angular_distance( - current_ra, current_dec, coverage_ras[ii + 1 :], coverage_decs[ii + 1 :] - ) - idx = np.where(dist <= 1 / 3600.0)[0] - if len(idx) > 0: - jj = idx[0] - diffs.append( - np.abs( - coverage_struct["data"][ii, 2] - coverage_struct["data"][jj, 2] - ) - ) - else: - for ii in range(len(coverage_struct["ipix"])): - ipix = coverage_struct["ipix"][ii] - for jj in range(len(coverage_struct["ipix"])): - if ii >= jj: - continue - if coverage_struct["telescope"][ii] == coverage_struct["telescope"][jj]: - continue - ipix2 = coverage_struct["ipix"][jj] - overlap = np.intersect1d(ipix, ipix2) - rat = np.array( - [ - float(len(overlap)) / float(len(ipix)), - float(len(overlap)) / float(len(ipix2)), - ] - ) - if np.any(rat > 0.5): - diffs.append( - np.abs( - coverage_struct["data"][ii, 2] - - coverage_struct["data"][jj, 2] - ) - ) - - export_tiles_coverage_hist(params, diffs) - - plot_tiles_coverage_hist(params, diffs) - plot_tiles_coverage_int( params, map_struct, catalog_struct, coverage_struct, plot_sun_moon ) diff --git a/gwemopt/plotting/plot_skymap.py b/gwemopt/plotting/plot_skymap.py index 14b56308..2a63592c 100644 --- a/gwemopt/plotting/plot_skymap.py +++ b/gwemopt/plotting/plot_skymap.py @@ -1,64 +1,30 @@ import healpy as hp +import ligo.skymap.plot import matplotlib.pyplot as plt import numpy as np +from astropy.io import fits +from ligo.skymap import moc +from matplotlib import pyplot as plt from gwemopt.plotting.style import CBAR_BOOL, UNIT, add_edges, cmap -def plot_skymap(params, map_struct): +def plot_skymap(params, map_struct, colnames=["PROB", "DISTMEAN", "DISTSTD"]): """ Function to plot the skymap """ plot_name = params["outputDir"].joinpath("prob.pdf") - if np.percentile(map_struct["prob"], 99) > 0: - hp.mollview( - map_struct["prob"], - title="", - unit=UNIT, - cbar=CBAR_BOOL, - min=np.percentile(map_struct["prob"], 1), - max=np.percentile(map_struct["prob"], 99), - cmap=cmap, - ) - else: - hp.mollview( - map_struct["prob"], - title="", - unit=UNIT, - cbar=CBAR_BOOL, - min=np.percentile(map_struct["prob"], 1), - cmap=cmap, - ) + hdu = map_struct["hdu"] + columns = [col.name for col in hdu.columns] - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() + for col in colnames: + if col in columns: - if "distmu" in map_struct: - fin = np.copy(map_struct["distmu"]) - fin[~np.isfinite(fin)] = np.nan - plot_name = params["outputDir"].joinpath("dist.pdf") - hp.mollview( - map_struct["distmu"], - unit="Distance [Mpc]", - min=np.nanpercentile(fin, 10), - max=np.nanpercentile(fin, 90), - ) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() - - fin = np.copy(map_struct["distmed"]) - fin[~np.isfinite(fin)] = np.nan - plot_name = params["outputDir"].joinpath("dist_median.pdf") - hp.mollview( - map_struct["distmed"], - unit="Distance [Mpc]", - min=np.nanpercentile(fin, 10), - max=np.nanpercentile(fin, 90), - ) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() + fig = plt.figure(figsize=(8, 6), dpi=100) + ax = plt.axes([0.05, 0.05, 0.9, 0.9], projection="astro mollweide") + ax.imshow_hpx(hdu, field=columns.index(col), cmap="cylon") + plot_name = params["outputDir"].joinpath(f"{col.lower()}.pdf") + plt.savefig(plot_name, dpi=200) + plt.close() diff --git a/gwemopt/plotting/plot_tiles.py b/gwemopt/plotting/plot_tiles.py index 7cc3736e..c0087893 100644 --- a/gwemopt/plotting/plot_tiles.py +++ b/gwemopt/plotting/plot_tiles.py @@ -3,6 +3,7 @@ import healpy as hp import matplotlib.pyplot as plt import numpy as np +from tqdm import tqdm from gwemopt.plotting.style import CBAR_BOOL, UNIT, add_edges, add_sun_moon, cmap @@ -13,62 +14,33 @@ def make_tile_plots(params, map_struct, tiles_structs, plot_sun_moon=True): """ plot_name = params["outputDir"].joinpath("tiles.pdf") - plt.figure() - hp.mollview(map_struct["prob"], title="", unit=UNIT, cbar=CBAR_BOOL, cmap=cmap) - ax = plt.gca() + + hdu = map_struct["hdu"] + columns = [col.name for col in hdu.columns] + + fig = plt.figure(figsize=(8, 6), dpi=100) + ax = plt.axes([0.05, 0.05, 0.9, 0.9], projection="astro mollweide") + ax.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon") for telescope in tiles_structs: tiles_struct = tiles_structs[telescope] - for index in tiles_struct.keys(): - patch = tiles_struct[index]["patch"] - if not patch: - continue - if type(patch) == list: - for p in patch: - hp.projaxes.HpxMollweideAxes.add_patch(ax, p) - else: - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch) + keys = list(tiles_struct.keys()) + probs = np.array([tiles_struct[key]["prob"] for key in keys]) + alphas = probs / np.max(probs) + + for ii, index in tqdm(enumerate(keys), total=len(keys)): + moc = tiles_struct[index]["moc"] + moc.fill( + ax=ax, + wcs=ax.wcs, + alpha=alphas[ii], + fill=True, + color="black", + linewidth=1, + ) + moc.border(ax=ax, wcs=ax.wcs, alpha=1, color="black") if plot_sun_moon: - add_sun_moon(params) + add_sun_moon(params, ax) - add_edges() plt.savefig(plot_name, dpi=200) plt.close() - - if params["doReferences"]: - rot = (90, 0, 0) - plot_name = params["outputDir"].joinpath("tiles_ref.pdf") - plt.figure() - hp.mollview( - np.zeros(map_struct["prob"].shape), - title="", - unit=UNIT, - cbar=CBAR_BOOL, - cmap=cmap, - ) - ax = plt.gca() - for telescope in tiles_structs: - config_struct = params["config"][telescope] - tiles_struct = tiles_structs[telescope] - if "reference_images" not in config_struct: - continue - for index in tiles_struct.keys(): - if index not in config_struct["reference_images"]: - continue - if len(params["filters"]) == 1: - if ( - not params["filters"][0] - in config_struct["reference_images"][index] - ): - continue - patch = tiles_struct[index]["patch"] - if not patch: - continue - patch_cpy = copy.copy(patch) - patch_cpy.axes = None - patch_cpy.figure = None - patch_cpy.set_transform(ax.transData) - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch_cpy) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() diff --git a/gwemopt/plotting/style.py b/gwemopt/plotting/style.py index 3275c146..d0a7788e 100644 --- a/gwemopt/plotting/style.py +++ b/gwemopt/plotting/style.py @@ -33,7 +33,7 @@ def add_edges(): hp.projtext(lon, lat, "%.0f" % lat, lonlat=True) -def add_sun_moon(params): +def add_sun_moon(params, ax): """ Add sun and moon position to the skymap """ @@ -44,11 +44,17 @@ def add_sun_moon(params): sun_position = get_body("sun", start_time) moon_position = get_body("moon", start_time) - hp.visufunc.projscatter( - sun_position.ra, sun_position.dec, lonlat=True, color="yellow" + ax.plot( + sun_position.ra, + sun_position.dec, + transform=ax.get_transform("world"), + color="yellow", ) - hp.visufunc.projscatter( - moon_position.ra, moon_position.dec, lonlat=True, color="grey" + ax.plot( + moon_position.ra, + moon_position.dec, + transform=ax.get_transform("world"), + color="grey", ) # also plot (in smaller scale) sun and moon position for the 7 following days @@ -62,23 +68,23 @@ def add_sun_moon(params): new_time = Time(new_time, format="gps", scale="utc") new_moon_position = get_body("moon", new_time) - hp.visufunc.projscatter( + ax.plot( new_moon_position.ra, new_moon_position.dec, - lonlat=True, + transform=ax.get_transform("world"), color="black", marker=".", - s=1, + markersize=1, ) if not i % 8: # only plot point for the sun every two days in order to avoid overlap new_sun_position = get_body("sun", new_time) - hp.visufunc.projscatter( + ax.plot( new_sun_position.ra, new_sun_position.dec, - lonlat=True, + transform=ax.get_transform("world"), color="black", marker=".", - s=1, + markersize=1, ) diff --git a/gwemopt/scheduler.py b/gwemopt/scheduler.py index ab4f8b29..26e52d01 100644 --- a/gwemopt/scheduler.py +++ b/gwemopt/scheduler.py @@ -438,9 +438,7 @@ def scheduler(params, config_struct, tile_struct): coverage_struct = {} coverage_struct["data"] = np.empty((0, 8)) coverage_struct["filters"] = [] - coverage_struct["ipix"] = [] - coverage_struct["patch"] = [] - coverage_struct["area"] = [] + coverage_struct["moc"] = [] if params["tilesType"] == "galaxy": coverage_struct["galaxies"] = [] @@ -541,13 +539,10 @@ def scheduler(params, config_struct, tile_struct): ) coverage_struct["filters"].append(filt) - coverage_struct["patch"].append(tile_struct_hold["patch"]) - coverage_struct["ipix"].append(tile_struct_hold["ipix"]) - coverage_struct["area"].append(tile_struct_hold["area"]) + coverage_struct["moc"].append(tile_struct_hold["moc"]) if params["tilesType"] == "galaxy": coverage_struct["galaxies"].append(tile_struct_hold["galaxies"]) - coverage_struct["area"] = np.array(coverage_struct["area"]) coverage_struct["filters"] = np.array(coverage_struct["filters"]) coverage_struct["FOV"] = [config_struct["FOV"]] * len(coverage_struct["filters"]) coverage_struct["telescope"] = [config_struct["telescope"]] * len( @@ -570,166 +565,3 @@ def computeSlewReadoutTime(config_struct, coverage_struct): prev_dec = dat[0] prev_ra = dat[1] return acc_time - - -def schedule_ra_splits( - params, - config_struct, - map_struct_hold, - tile_struct, - telescope, - previous_coverage_struct, -): - location = astropy.coordinates.EarthLocation( - config_struct["longitude"], - config_struct["latitude"], - config_struct["elevation"], - ) - - raslices = gwemopt.utils.utils.auto_rasplit( - params, map_struct_hold, params["nside_down"] - ) - - maxidx = 0 - coverage_structs = [] - skip = False - while len(raslices) != 0: - params["unbalanced_tiles"] = [] - config_struct["exposurelist"] = segments.segmentlist( - config_struct["exposurelist"][maxidx:] - ) - if len(config_struct["exposurelist"]) < 2: - break - - map_struct_slice = copy.deepcopy(map_struct_hold) - - exposurelist = np.array_split(config_struct["exposurelist"], len(raslices))[0] - minhas = [] - minhas_late = [] - try_end = False - if len(raslices) == 1: - raslice = raslices[0] - del raslices[0] - else: - for raslice in raslices: - has = [] - has_late = [] - for seg in exposurelist: - mjds = np.linspace(seg[0], seg[1], 100) - tt = Time(mjds, format="mjd", scale="utc", location=location) - lst = tt.sidereal_time("mean") / u.hourangle - ha = np.abs(lst - raslice[0]) - ha_late = np.abs(lst - raslice[1]) - - idx = np.where(ha > 12.0)[0] - ha[idx] = 24.0 - ha[idx] - idx_late = np.where(ha_late > 12.0)[0] - ha_late[idx_late] = 24.0 - ha_late[idx_late] - has += list(ha) - has_late += list(ha_late) - if len(has) > 0: - minhas.append(np.min(has)) - if len(has_late) > 0: - minhas_late.append(np.min(has_late)) - - if (len(minhas_late) > 0) and (len(has_late) > 0): - # conditions for trying to schedule end of slice - if np.min(minhas_late) <= 5.0 and np.min(has) > 4.0 and not skip: - try_end = True - min = np.argmin(minhas_late) - raslice = raslices[min] - else: - min = np.argmin(minhas) - raslice = raslices[min] - del raslices[min] - else: - min = np.argmin(minhas) - raslice = raslices[min] - del raslices[min] - - # do RA slicing - ra_low, ra_high = raslice[0], raslice[1] - ra = map_struct_slice["ra"] - if ra_low <= ra_high: - ipix = np.where( - (ra_high * 360.0 / 24.0 < ra) | (ra_low * 360.0 / 24.0 > ra) - )[0] - else: - ipix = np.where( - (ra_high * 360.0 / 24.0 < ra) & (ra_low * 360.0 / 24.0 > ra) - )[0] - - map_struct_slice["prob"][ipix] = 0.0 - - if params["timeallocationType"] == "absmag": - tile_struct = gwemopt.tiles.absmag_tiles_struct( - params, config_struct, telescope, map_struct_slice, tile_struct - ) - else: - tile_struct = gwemopt.tiles.powerlaw_tiles_struct( - params, config_struct, telescope, map_struct_slice, tile_struct - ) - - config_struct_hold = copy.copy(config_struct) - coverage_struct, tile_struct = schedule_alternating( - params, - config_struct_hold, - telescope, - map_struct_slice, - tile_struct, - previous_coverage_struct, - ) - if len(coverage_struct["ipix"]) == 0: - continue - optimized_max, coverage_struct, tile_struct = optimize_max_tiles( - params, - tile_struct, - coverage_struct, - config_struct, - telescope, - map_struct_slice, - ) - params["max_nb_tiles"] = np.array([optimized_max], dtype=float) - balanced_fields = 0 - coverage_struct, tile_struct = schedule_alternating( - params, - config_struct, - telescope, - map_struct_slice, - tile_struct, - previous_coverage_struct, - ) - - doReschedule, balanced_fields = balance_tiles( - params, tile_struct, coverage_struct - ) - config_struct_hold = copy.copy(config_struct) - - if balanced_fields == 0: - if try_end: - skip = True - continue - elif try_end: - del raslices[min] - skip = False - - if len(coverage_struct["exposureused"]) > 0: - maxidx = int(coverage_struct["exposureused"][-1]) - - coverage_struct = gwemopt.utils.utils.erase_unbalanced_tiles( - params, coverage_struct - ) - - # limit to max number of filter sets - if len(coverage_structs) < params["max_filter_sets"]: - coverage_structs.append(coverage_struct) - else: - prob_structs = [ - np.sum(prev_struct["data"][:, 6]) for prev_struct in coverage_structs - ] - if np.any(np.array(prob_structs) < np.sum(coverage_struct["data"][:, 6])): - argmin = np.argmin(prob_structs) - del coverage_structs[argmin] - coverage_structs.append(coverage_struct) - - return gwemopt.coverage.combine_coverage_structs(coverage_structs) diff --git a/gwemopt/tiles.py b/gwemopt/tiles.py index a9aaa7cb..056b261e 100644 --- a/gwemopt/tiles.py +++ b/gwemopt/tiles.py @@ -1,12 +1,15 @@ +import bisect import copy +import astropy_healpix as ah import healpy as hp import ligo.segments as segments import numpy as np import pandas as pd from astropy import units as u -from astropy.coordinates import SkyCoord +from astropy.coordinates import ICRS, SkyCoord from astropy.time import Time +from mocpy import MOC from scipy.stats import norm from shapely.geometry import MultiPoint from tqdm import tqdm @@ -16,6 +19,10 @@ import gwemopt.segments from gwemopt.utils.geometry import angular_distance +LEVEL = MOC.MAX_ORDER +HPX = ah.HEALPix(nside=ah.level_to_nside(LEVEL), order="nested", frame=ICRS()) +PIXEL_AREA = HPX.pixel_area.to_value(u.sr) + TILE_TYPES = ["moc", "galaxy"] @@ -1071,9 +1078,9 @@ def powerlaw_tiles_struct( tot_obs_time = config_struct["tot_obs_time"] if "observability" in map_struct: - prob = map_struct["observability"][telescope]["prob"] + skymap = map_struct["observability"][telescope]["prob"] else: - prob = map_struct["prob"] + skymap = map_struct["skymap"] n, cl, dist_exp = ( params["powerlaw_n"], @@ -1085,96 +1092,29 @@ def powerlaw_tiles_struct( tile_probs = compute_tiles_map( params, tile_struct, - prob, + skymap, func="galaxy", - ipix_keep=map_struct["ipix_keep"], + moc_keep=map_struct["moc_keep"], catalog_struct=catalog_struct, ) else: tile_probs = compute_tiles_map( params, tile_struct, - prob, + skymap, func="np.sum(x)", - ipix_keep=map_struct["ipix_keep"], + moc_keep=map_struct["moc_keep"], ) tile_probs[tile_probs < np.max(tile_probs) * 0.01] = 0.0 - prob_scaled = copy.deepcopy(prob) - prob_sorted = np.sort(prob_scaled)[::-1] - prob_indexes = np.argsort(prob_scaled)[::-1] - prob_cumsum = np.cumsum(prob_sorted) - index = np.argmin(np.abs(prob_cumsum - cl)) + 1 - # prob_scaled[prob_indexes[index:]] = 1e-10 - prob_scaled[prob_indexes[index:]] = 0.0 - prob_scaled = prob_scaled**n - prob_scaled = prob_scaled / np.nansum(prob_scaled) - - if params["tilesType"] == "galaxy": - ranked_tile_probs = compute_tiles_map( - params, - tile_struct, - prob_scaled, - func="galaxy", - ipix_keep=map_struct["ipix_keep"], - catalog_struct=catalog_struct, - ) - else: - ranked_tile_probs = compute_tiles_map( - params, - tile_struct, - prob_scaled, - func="np.sum(x)", - ipix_keep=map_struct["ipix_keep"], - ) - - ranked_tile_probs[np.isnan(ranked_tile_probs)] = 0.0 - ranked_tile_probs_thresh = np.max(ranked_tile_probs) * 0.01 - ranked_tile_probs[ranked_tile_probs <= ranked_tile_probs_thresh] = 0.0 - ranked_tile_probs = ranked_tile_probs / np.nansum(ranked_tile_probs) - - if "distmed" in map_struct: - distmed = map_struct["distmed"] - distmed[distmed <= 0] = np.nan - distmed[~np.isfinite(distmed)] = np.nan - # distmed[distmed 0, ii] = params["exposuretimes"][ii] + ranked_tile_times[tile_probs > 0, ii] = params["exposuretimes"][ii] - for key, prob, exposureTime, tileprob in zip( - keys, ranked_tile_probs, ranked_tile_times, tile_probs - ): + for key, tileprob, exposureTime in zip(keys, tile_probs, ranked_tile_times): # Try to load the minimum duration of time from telescope config file # Otherwise set it to zero try: @@ -1203,7 +1143,7 @@ def powerlaw_tiles_struct( tileprob = 0.0 tile_struct[key]["prob"] = tileprob - if prob == 0.0: + if tileprob == 0.0: tile_struct[key]["exposureTime"] = 0.0 tile_struct[key]["nexposures"] = 0 tile_struct[key]["filt"] = [] @@ -1236,7 +1176,7 @@ def powerlaw_tiles_struct( else: ranked_tile_times = gwemopt.utils.integrationTime( tot_obs_time, - ranked_tile_probs, + tile_probs, func=None, T_int=config_struct["exposuretime"], ) @@ -1244,7 +1184,7 @@ def powerlaw_tiles_struct( keys = tile_struct.keys() for key, prob, exposureTime, tileprob in zip( - keys, ranked_tile_probs, ranked_tile_times, tile_probs + keys, tile_probs, ranked_tile_times ): # Try to load the minimum duration of time from telescope config file # Otherwise set it to zero @@ -1309,90 +1249,44 @@ def moc(params, map_struct, moc_structs, doSegments=True): return tile_structs -def pem_tiles_struct(params, config_struct, telescope, map_struct, tile_struct): - tot_obs_time = config_struct["tot_obs_time"] - - if "observability" in map_struct: - prob = map_struct["observability"][telescope]["prob"] - else: - prob = map_struct["prob"] - - ranked_tile_probs = compute_tiles_map( - params, tile_struct, prob, func="np.sum(x)", ipix_keep=map_struct["ipix_keep"] - ) - - lim_mag = config_struct["magnitude"] - lim_time = config_struct["exposuretime"] +def find_overlap(intervals_list, query_intervals, probdensity): - tau = np.arange(lim_time, 3600.0, lim_time) - Loftau = None + # Sort both intervals list and query_intervals list based on start points + intervals_list.sort(key=lambda x: x[0]) + t0s = np.array([query_interval[0] for query_interval in query_intervals]) + idx = np.argsort(t0s).tolist() + query_intervals = [query_intervals[i] for i in idx] + probdensity = [probdensity[i] for i in idx] - N_ref = 9.7847e9 - L_min = 4.9370e31 - L_max = 4.9370e33 - - if "distmu" in map_struct: - prob = map_struct["prob"] - distnorm = map_struct["distnorm"] - distmu = map_struct["distmu"] - distsigma = map_struct["distsigma"] - - D_min = 1.0e7 - D_max = 1.0e10 - R = np.linspace(D_min / 1e6, D_max / 1e6) - p_R = [ - np.sum(prob * rr**2 * distnorm * norm(distmu, distsigma).pdf(rr)) - for rr in R - ] - p_R = np.array(p_R) - p_R = p_R / len(prob) - - R = R * 1e6 - D_mu = None - D_sig = None - else: - D_mu = 200.0e6 - D_sig = 60.0e6 - R = None - p_R = None - - tau, prob = gwemopt.pem.Pem( - lim_mag, - lim_time, - N_ref=N_ref, - L_min=L_min, - L_max=L_max, - tau=tau, - Loftau=Loftau, - D_mu=D_mu, - D_sig=D_sig, - R=R, - p_R=p_R, - ) - - tprob, time_allocation = gwemopt.pem.Main( - tot_obs_time, 0, 0, ranked_tile_probs, tau, prob, "Eq" - ) - - if params["doPlots"]: - gwemopt.plotting.tauprob(params, tau, prob) - - keys = tile_struct.keys() - for key, prob, exposureTime in zip(keys, ranked_tile_probs, time_allocation): - tile_struct[key]["prob"] = prob - tile_struct[key]["exposureTime"] = exposureTime - tile_struct[key]["nexposures"] = int( - np.floor(exposureTime / config_struct["exposuretime"]) - ) - tile_struct[key]["filt"] = [config_struct["filt"]] * tile_struct[key][ - "nexposures" - ] + prob = 0 + for interval in intervals_list: + interval_start, interval_end = interval + # Find the index where interval_start would be inserted in query_intervals + index = bisect.bisect_right(query_intervals, [interval_start, float("inf")]) + # Iterate over query_intervals starting from the index to find overlaps + for query_interval, probdens in zip( + query_intervals[index:], probdensity[index:] + ): + query_start, query_end = query_interval + if query_start > interval_end: + break # No more overlaps possible, break the loop + overlap_start = max(interval_start, query_start) + overlap_end = min(interval_end, query_end) + overlap = max(0, overlap_end - overlap_start) + if overlap > 0: + # length * area * overlap + prob = prob + overlap * PIXEL_AREA * probdens - return tile_struct + return prob def compute_tiles_map( - params, tile_struct, skymap, func=None, ipix_keep=[], catalog_struct=None + params, + tile_struct, + skymap, + func=None, + moc_keep=MOC.new_empty(max_depth=10), + catalog_struct=None, ): if func is None: f = lambda x: np.sum(x) @@ -1426,31 +1320,12 @@ def compute_tiles_map( keys = tile_struct.keys() ntiles = len(keys) vals = np.nan * np.ones((ntiles,)) - for ii, key in enumerate(tile_struct.keys()): - idx = np.where(prob[tile_struct[key]["ipix"]] == -1)[0] - idx = np.setdiff1d(idx, ipix_keep) - if len(prob[tile_struct[key]["ipix"]]) == 0: - rat = 0.0 - else: - rat = float(len(idx)) / float(len(prob[tile_struct[key]["ipix"]])) - if rat > params["maximumOverlap"]: - vals[ii] = 0.0 - else: - ipix = tile_struct[key]["ipix"] - if len(ipix) == 0: - vals[ii] = 0.0 - else: - vals_to_sum = prob[ipix] - if func == "np.nanmedian(x)": - vals_to_sum[vals_to_sum < 0] = np.nan - else: - vals_to_sum[vals_to_sum < 0] = 0 - vals[ii] = f(vals_to_sum) - - ipix_keep = np.setdiff1d(ipix_keep, tile_struct[key]["ipix"]) - ipix_slice = np.setdiff1d(tile_struct[key]["ipix"], ipix_keep) - if len(ipix_slice) > 0: - prob[ipix_slice] = -1 + for ii, key in tqdm(enumerate(keys), total=len(keys)): + mr1 = [] + for lo, hi in tile_struct[key]["moc"].to_depth29_ranges: + mr1.append([lo, hi]) + mr2 = skymap["HPX"].data.tolist() + vals[ii] = find_overlap(mr1, mr2, skymap["PROBDENSITY"].tolist()) return vals diff --git a/gwemopt/utils/pixels.py b/gwemopt/utils/pixels.py index 415a7eb9..6b14b684 100644 --- a/gwemopt/utils/pixels.py +++ b/gwemopt/utils/pixels.py @@ -213,35 +213,21 @@ def getRectanglePixels( idx1 = np.where(np.abs(radecs[:, 1]) >= 87.0)[0] if len(idx1) == 4: - return [], [], [], [] + return MOC.new_empty(max_depth=10) idx1 = np.where((radecs[:, 1] >= 87.0) | (radecs[:, 1] <= -87.0))[0] if len(idx1) > 0: radecs = np.delete(radecs, idx1[0], 0) - xyz = [] - for r, d in radecs: - xyz.append(hp.ang2vec(r, d, lonlat=True)) - npts, junk = radecs.shape if npts == 4: - xyz = [xyz[0], xyz[1], xyz[3], xyz[2]] - ipix = hp.query_polygon(nside, np.array(xyz)) - else: - ipix = hp.query_polygon(nside, np.array(xyz)) + idx = [0, 1, 3, 2] + radecs = radecs[idx, :] - xyz = np.array(xyz) - proj = hp.projector.MollweideProj(rot=rotation, coord=None) - x, y = proj.vec2xy(xyz[:, 0], xyz[:, 1], xyz[:, 2]) - xy = np.zeros(radecs.shape) - xy[:, 0] = x - xy[:, 1] = y - path = matplotlib.path.Path(xy) - patch = matplotlib.patches.PathPatch( - path, alpha=alpha, facecolor=color, fill=True, zorder=3, edgecolor=edgecolor - ) + coords = coordinates.SkyCoord(radecs[:, 0] * u.deg, radecs[:, 1] * u.deg) + moc = MOC.from_polygon_skycoord(coords) - return ipix, radecs, patch, area + return moc def getSquarePixels( @@ -289,32 +275,18 @@ def getSquarePixels( idx1 = np.where(np.abs(radecs[:, 1]) >= 87.0)[0] if len(idx1) == 4: - return [], [], [], [] + return MOC.new_empty(max_depth=10) idx1 = np.where((radecs[:, 1] >= 87.0) | (radecs[:, 1] <= -87.0))[0] if len(idx1) > 0: radecs = np.delete(radecs, idx1[0], 0) - xyz = [] - for r, d in radecs: - xyz.append(hp.ang2vec(r, d, lonlat=True)) - npts, junk = radecs.shape if npts == 4: - xyz = [xyz[0], xyz[1], xyz[3], xyz[2]] - ipix = hp.query_polygon(nside, np.array(xyz)) - else: - ipix = hp.query_polygon(nside, np.array(xyz)) + idx = [0, 1, 3, 2] + radecs = radecs[idx, :] - xyz = np.array(xyz) - proj = hp.projector.MollweideProj(rot=rotation, coord=None) - x, y = proj.vec2xy(xyz[:, 0], xyz[:, 1], xyz[:, 2]) - xy = np.zeros(radecs.shape) - xy[:, 0] = x - xy[:, 1] = y - path = matplotlib.path.Path(xy) - patch = matplotlib.patches.PathPatch( - path, alpha=alpha, facecolor=color, fill=True, zorder=3, edgecolor=edgecolor - ) + coords = coordinates.SkyCoord(radecs[:, 0] * u.deg, radecs[:, 1] * u.deg) + moc = MOC.from_polygon_skycoord(coords) - return ipix, radecs, patch, area + return moc