diff --git a/.github/workflows/black.yml b/.github/workflows/black.yml index ab3e0d0..fcf59af 100644 --- a/.github/workflows/black.yml +++ b/.github/workflows/black.yml @@ -14,7 +14,7 @@ jobs: - uses: actions/setup-python@v4 with: - python-version: ${{ matrix.python-version }} + python-version: '3.12' cache: 'pip' cache-dependency-path: setup.py diff --git a/.github/workflows/continous_integration.yml b/.github/workflows/continous_integration.yml index 55a780a..a865b9d 100644 --- a/.github/workflows/continous_integration.yml +++ b/.github/workflows/continous_integration.yml @@ -32,24 +32,34 @@ jobs: steps: # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - - uses: actions/setup-python@v4 + - uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} cache: 'pip' cache-dependency-path: setup.py + - uses: actions/cache@v4 + with: + path: | + ~/.cache/pip + key: ${{ runner.os }}-pip-${{ hashFiles('**/pyproject.toml') }} + restore-keys: | + ${{ runner.os }}-pip-${{ hashFiles('**/pyproject.toml') }} + ${{ runner.os }}-pip- + - name: Install dependencies run: | sudo apt-get update sudo apt-get install -y ffmpeg python -m pip install --upgrade pip - pip install pytest pytest-runner pytest-rerunfailures coveralls freezegun - pip install -e . + pip install -e ".[test]" - name: Run tests - run: coverage run -m pytest -v -s + run: | + export NUMBA_DISABLE_JIT=1 + python -m coverage run -m pytest -v -s - name: Run Coveralls if: ${{ success() }} diff --git a/.github/workflows/isort.yml b/.github/workflows/isort.yml index ffd3206..d7e5aa2 100644 --- a/.github/workflows/isort.yml +++ b/.github/workflows/isort.yml @@ -9,10 +9,10 @@ jobs: isort: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 - - uses: actions/setup-python@v4 + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 with: - python-version: ${{ matrix.python-version }} + python-version: '3.12' cache: 'pip' cache-dependency-path: setup.py diff --git a/README.md b/README.md index b469b91..bc5f950 100644 --- a/README.md +++ b/README.md @@ -62,11 +62,12 @@ If you want the latest version, we recommend creating a clean environment: ```commandline conda create -n gwemopt python=3.11 git clone git@github.com:skyportal/gwemopt.git -pip install -e gwemopt +cd gwemopt +pip install -e . pre-commit install ``` -or if you just want the latest version on Github: +or if you just want the latest version released on PyPI: ```commandline pip install gwemopt @@ -80,7 +81,7 @@ conda install -c astropy astroquery conda install -c conda-forge voeventlib astropy-healpix python-ligo-lw ligo-segments ligo.skymap ffmpeg ``` -And then run `pip install -e gwemopt` again. +And then run `pip install -e .` again. # Usage @@ -91,3 +92,17 @@ gwemopt-run .... ``` where ... corresponds to the various arguments. + +# Tests + +To run the tests, you'll first need to install gwemopt with the testing dependencies: + +```commandline +pip install -e ".[test]" +``` + +Then you can run the tests using: + +```commandline +coverage run -m pytest -v -s +``` diff --git a/gwemopt/catalogs/get.py b/gwemopt/catalogs/get.py index 607a94c..f63647a 100644 --- a/gwemopt/catalogs/get.py +++ b/gwemopt/catalogs/get.py @@ -2,9 +2,7 @@ import healpy as hp import numpy as np -from astropy.table import Table from astroquery.vizier import Vizier -from ligo.skymap import distance from ligo.skymap.bayestar import derasterize from scipy.stats import norm diff --git a/gwemopt/catalogs/mangrove.py b/gwemopt/catalogs/mangrove.py index b0439bb..c37ae13 100644 --- a/gwemopt/catalogs/mangrove.py +++ b/gwemopt/catalogs/mangrove.py @@ -1,5 +1,4 @@ import subprocess -from pathlib import Path import numpy as np import pandas as pd diff --git a/gwemopt/catalogs/nedlvs.py b/gwemopt/catalogs/nedlvs.py index 434ce4c..8750483 100644 --- a/gwemopt/catalogs/nedlvs.py +++ b/gwemopt/catalogs/nedlvs.py @@ -1,5 +1,4 @@ import subprocess -from pathlib import Path import numpy as np import pandas as pd diff --git a/gwemopt/catalogs/twomrs.py b/gwemopt/catalogs/twomrs.py index 8fb7ecc..37de70a 100644 --- a/gwemopt/catalogs/twomrs.py +++ b/gwemopt/catalogs/twomrs.py @@ -4,6 +4,7 @@ from astropy import units as u from astropy.cosmology import WMAP9 as cosmo from astropy.io.misc.hdf5 import read_table_hdf5, write_table_hdf5 +from astropy.table import Table from astroquery.vizier import Vizier from scipy.special import gammaincinv diff --git a/gwemopt/chipgaps/decam.py b/gwemopt/chipgaps/decam.py index 20e8e6b..7043ee9 100644 --- a/gwemopt/chipgaps/decam.py +++ b/gwemopt/chipgaps/decam.py @@ -18,15 +18,12 @@ This module contains the DECamtile class, which is used to calculate decam chip gaps """ -from array import array - import astropy -import healpy as hp import numpy as np from astropy import units as u from astropy.coordinates import ICRS, SkyCoord from astropy.wcs import WCS -from mocpy import MOC, mocpy +from mocpy import MOC from tqdm import tqdm @@ -176,7 +173,6 @@ def get_decam_ccds(ra, dec, save_footprint=False): # ccd_prob = CCDProb(ra, dec) decam_tile = DECamtile(ra, dec) ccd_cents_ra = decam_tile.ccd_RA - ccd_cents_dec = decam_tile.ccd_Dec offsets = np.asarray( [ decam_tile._wcs[ccd_id].calc_footprint(axes=decam_tile.ccd_size) @@ -229,10 +225,7 @@ def get_decam_quadrant_moc(ra, dec, n_threads=None): stacked = np.stack((ccd_coords.ra.deg, ccd_coords.dec.deg), axis=1) result = stacked.reshape(-1, ccd_coords.ra.deg.shape[1]) lon_lat_list = [row for row in result] - moc = sum( - MOC.from_polygons(lon_lat_list, max_depth=np.uint8(10), n_threads=n_threads) - ) - mocs.append(moc) + mocs.append(sum(MOC.from_polygons(lon_lat_list, False, 10, n_threads))) return mocs diff --git a/gwemopt/chipgaps/ztf.py b/gwemopt/chipgaps/ztf.py index 9485715..52f8248 100644 --- a/gwemopt/chipgaps/ztf.py +++ b/gwemopt/chipgaps/ztf.py @@ -20,7 +20,6 @@ """ import astropy -import healpy as hp import numpy as np from astropy import units as u from astropy.coordinates import ICRS, SkyCoord @@ -290,9 +289,6 @@ def get_ztf_quadrant_moc(ra, dec, n_threads=None): stacked = np.stack((ccd_coords.ra.deg, ccd_coords.dec.deg), axis=1) result = stacked.reshape(-1, ccd_coords.ra.deg.shape[1]) lon_lat_list = [row for row in result] - moc = sum( - MOC.from_polygons(lon_lat_list, max_depth=np.uint8(10), n_threads=n_threads) - ) - mocs.append(moc) + mocs.append(sum(MOC.from_polygons(lon_lat_list, False, 10, n_threads))) return mocs diff --git a/gwemopt/coverage.py b/gwemopt/coverage.py index 9ef925d..c6a052a 100644 --- a/gwemopt/coverage.py +++ b/gwemopt/coverage.py @@ -1,13 +1,11 @@ import copy -import ligo.segments as segments import numpy as np import regions from astropy import units as u from astropy.coordinates import SkyCoord, get_sun from astropy.time import Time from mocpy import MOC -from tqdm import tqdm import gwemopt.plotting import gwemopt.scheduler @@ -16,7 +14,6 @@ from gwemopt.tiles import ( balance_tiles, check_overlapping_tiles, - eject_tiles, optimize_max_tiles, order_by_observability, slice_galaxy_tiles, @@ -60,7 +57,6 @@ def combine_coverage_structs(coverage_structs): def read_coverage(params, telescope, filename, moc_struct=None): - nside = params["nside"] config_struct = params["config"][telescope] schedule_table = read_schedule(filename) @@ -69,7 +65,7 @@ def read_coverage(params, telescope, filename, moc_struct=None): coverage_struct["filters"] = [] coverage_struct["moc"] = [] - for ii, row in schedule_table.iterrows(): + for _, row in schedule_table.iterrows(): ra, dec = row["ra"], row["dec"] tobs = row["tobs"] limmag = row["limmag"] @@ -139,10 +135,9 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): map_struct_hold = copy.deepcopy(map_struct) coverage_structs = [] - n_scope = 0 full_prob_map = map_struct["skymap"] - for jj, telescope in enumerate(params["telescopes"]): + for telescope in params["telescopes"]: config_struct = params["config"][telescope] tile_struct = tile_structs[telescope] @@ -167,7 +162,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): params_hold = copy.copy(params) config_struct_hold = copy.copy(config_struct) - coverage_struct, tile_struct = gwemopt.scheduler.schedule_alternating( + coverage_struct, tile_struct = gwemopt.tiles.schedule_alternating( params_hold, config_struct_hold, telescope, @@ -196,7 +191,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): ( coverage_struct, tile_struct, - ) = gwemopt.scheduler.schedule_alternating( + ) = gwemopt.tiles.schedule_alternating( params_hold, config_struct_hold, telescope, @@ -204,16 +199,14 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): tile_struct, previous_coverage_struct, ) - doReschedule, balanced_fields = balance_tiles( - params_hold, tile_struct, coverage_struct - ) + doReschedule, _ = balance_tiles(params_hold, coverage_struct) if doReschedule: config_struct_hold = copy.copy(config_struct) ( coverage_struct, tile_struct, - ) = gwemopt.scheduler.schedule_alternating( + ) = gwemopt.tiles.schedule_alternating( params_hold, config_struct_hold, telescope, @@ -334,9 +327,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): if params["doBalanceExposure"]: cnt, ntrials = 0, 10 while cnt < ntrials: - doReschedule, balanced_fields = balance_tiles( - params, tile_struct, coverage_struct - ) + doReschedule, _ = balance_tiles(params, coverage_struct) if doReschedule: for key in params["unbalanced_tiles"]: tile_struct[key]["prob"] = 0.0 @@ -351,7 +342,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): if params["max_nb_tiles"] is not None: tile_struct, doReschedule = slice_number_tiles( - params, telescope, tile_struct, coverage_struct + params, tile_struct, coverage_struct ) if doReschedule: coverage_struct = gwemopt.scheduler.scheduler( diff --git a/gwemopt/efficiency.py b/gwemopt/efficiency.py index 112cdce..e484ca7 100644 --- a/gwemopt/efficiency.py +++ b/gwemopt/efficiency.py @@ -6,7 +6,6 @@ import scipy.stats from astropy.coordinates import SkyCoord from astropy.time import Time -from ligo.skymap import distance from scipy.interpolate import interp1d from gwemopt.io.export_efficiency import ( @@ -94,9 +93,7 @@ def compute_efficiency(params, map_struct, lightcurve_struct, coverage_struct): efficiency_struct["dists_inj"] = dists_inj if params["true_location"]: - det = compute_true_efficiency( - params, map_struct, lightcurve_struct, coverage_struct - ) + det = compute_true_efficiency(params, lightcurve_struct, coverage_struct) lc_name = lightcurve_struct["name"] save_detection_nondetection( @@ -108,10 +105,7 @@ def compute_efficiency(params, map_struct, lightcurve_struct, coverage_struct): return efficiency_struct -def compute_true_efficiency(params, map_struct, lightcurve_struct, coverage_struct): - nside = params["nside"] - npix = hp.nside2npix(nside) - Ninj = params["Ninj"] +def compute_true_efficiency(params, lightcurve_struct, coverage_struct): gpstime = params["gpstime"] mjd_inj = Time(gpstime, format="gps", scale="utc").mjd @@ -133,9 +127,7 @@ def compute_true_efficiency(params, map_struct, lightcurve_struct, coverage_stru lightcurve_mag = lightcurve_struct[filt] idx = np.where(np.isfinite(lightcurve_mag))[0] - f = interp.interp1d( - lightcurve_t[idx], lightcurve_mag[idx], fill_value="extrapolate" - ) + f = interp1d(lightcurve_t[idx], lightcurve_mag[idx], fill_value="extrapolate") lightcurve_mag_interp = f(mjd) dist_threshold = (10 ** (((mag - lightcurve_mag_interp) / 5.0) + 1.0)) / 1e6 diff --git a/gwemopt/io/schedule.py b/gwemopt/io/schedule.py index 3f5365e..fa3d85a 100644 --- a/gwemopt/io/schedule.py +++ b/gwemopt/io/schedule.py @@ -60,7 +60,7 @@ def read_schedule(schedule_path): return schedule -def export_schedule_xml(xmlfile, map_struct, coverage_struct, config_struct): +def export_schedule_xml(xmlfile, coverage_struct, config_struct): what = What() table = Table(name="data", Description=["The datas of GWAlert"]) @@ -157,9 +157,8 @@ def export_schedule_xml(xmlfile, map_struct, coverage_struct, config_struct): data = coverage_struct["data"][ii, :] ra, dec = data[0], data[1] - observ_time, exposure_time, field_id, prob, airmass = ( + observ_time, field_id, prob, airmass = ( data[2], - data[4], data[5], data[6], data[7], @@ -195,7 +194,6 @@ def export_schedule_xml(xmlfile, map_struct, coverage_struct, config_struct): def summary(params, map_struct, coverage_struct, catalog_struct=None): - filts = list(set(coverage_struct["filters"])) for jj, telescope in enumerate(params["telescopes"]): schedulefile = params["outputDir"].joinpath(f"schedule_{telescope}.dat") @@ -203,7 +201,7 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): config_struct = params["config"][telescope] - export_schedule_xml(schedulexmlfile, map_struct, coverage_struct, config_struct) + export_schedule_xml(schedulexmlfile, coverage_struct, config_struct) if (params["tilesType"] == "hierarchical") or (params["tilesType"] == "greedy"): fields = np.zeros((params["Ntiles"][jj], len(filts) + 2)) @@ -294,7 +292,6 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): for tt in tts: mjds_floor = [] mjds = [] - ipixs = np.empty((0, 2)) cum_prob = 0.0 cum_area = 0.0 diff --git a/gwemopt/io/skymap.py b/gwemopt/io/skymap.py index 7762b6d..bb39444 100644 --- a/gwemopt/io/skymap.py +++ b/gwemopt/io/skymap.py @@ -8,7 +8,6 @@ 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 @@ -20,8 +19,6 @@ from ligo.skymap import moc from ligo.skymap.bayestar import derasterize, rasterize from ligo.skymap.io import read_sky_map -from matplotlib import pyplot as plt -from mocpy import MOC from scipy.interpolate import PchipInterpolator from scipy.stats import norm @@ -291,11 +288,7 @@ def read_skymap(params, map_struct=None): 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 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 @@ -334,7 +327,7 @@ def read_skymap(params, map_struct=None): ( map_struct["skymap_raster"]["DISTMEAN"], map_struct["skymap_raster"]["DISTSTD"], - mom_norm, + _, ) = ligodist.parameters_to_moments( map_struct["skymap_raster"]["DISTMU"], map_struct["skymap_raster"]["DISTSIGMA"], diff --git a/gwemopt/moc.py b/gwemopt/moc.py index 56b9e5d..7892dcc 100644 --- a/gwemopt/moc.py +++ b/gwemopt/moc.py @@ -2,23 +2,17 @@ import os import astropy.units as u -import healpy as hp import numpy as np import regions -from astropy.coordinates import SkyCoord -from joblib import Parallel, delayed from mocpy import MOC -from tqdm import tqdm import gwemopt.tiles from gwemopt.chipgaps import get_decam_quadrant_moc, get_ztf_quadrant_moc from gwemopt.paths import CONFIG_DIR -from gwemopt.utils.parallel import tqdm_joblib from gwemopt.utils.pixels import get_region_moc def construct_moc(params, config_struct, telescope, tesselation): - if params["doParallel"]: n_threads = params["Ncores"] else: @@ -42,7 +36,7 @@ def construct_moc(params, config_struct, telescope, tesselation): lon=tesselation[:, 1] * u.deg, lat=tesselation[:, 2] * u.deg, radius=config_struct["FOV"] * u.deg, - max_depth=np.uint8(10), + max_depth=10, n_threads=n_threads, ) elif config_struct["FOV_type"] == "square": @@ -52,7 +46,7 @@ def construct_moc(params, config_struct, telescope, tesselation): a=config_struct["FOV"] * u.deg, b=config_struct["FOV"] * u.deg, angle=0 * u.deg, - max_depth=np.uint8(10), + max_depth=10, n_threads=n_threads, ) elif config_struct["FOV_type"] == "region": @@ -69,8 +63,6 @@ def construct_moc(params, config_struct, telescope, tesselation): def create_moc(params, map_struct=None, field_ids=None, from_skyportal=False): - nside = params["nside"] - moc_structs = {} for telescope in params["telescopes"]: config_struct = params["config"][telescope] @@ -105,7 +97,6 @@ def create_moc(params, map_struct=None, field_ids=None, from_skyportal=False): moc_struct[index]["ra"] = tess.ra moc_struct[index]["dec"] = tess.dec moc_struct[index]["moc"] = moc - else: if (telescope == "ZTF") and params["doUsePrimary"]: idx = np.where(tesselation[:, 0] <= 880)[0] diff --git a/gwemopt/plotting/plot_coverage.py b/gwemopt/plotting/plot_coverage.py index a5f7ddc..6abbda0 100644 --- a/gwemopt/plotting/plot_coverage.py +++ b/gwemopt/plotting/plot_coverage.py @@ -1,5 +1,3 @@ -import copy - import healpy as hp import matplotlib.pyplot as plt import numpy as np @@ -10,7 +8,6 @@ from gwemopt.io import export_tiles_coverage_int from gwemopt.plotting.movie import make_movie from gwemopt.plotting.style import add_sun_moon -from gwemopt.utils.geometry import angular_distance def plot_tiles_coverage(params, map_struct, coverage_struct, plot_sun_moon=False): @@ -266,7 +263,6 @@ def plot_coverage_scaled(params, map_struct, coverage_struct, plot_sun_moon, max def plot_coverage_movie(params, map_struct, coverage_struct, plot_sun_moon, max_time): - hdu = map_struct["hdu"] columns = [col.name for col in hdu.columns] @@ -330,7 +326,6 @@ def make_plot(jj, params, map_struct, coverage_struct): def make_coverage_plots( params, map_struct, coverage_struct, catalog_struct=None, plot_sun_moon: bool = True ): - idx = np.isfinite(coverage_struct["data"][:, 4]) if not idx.size: return @@ -347,7 +342,6 @@ def make_coverage_plots( def make_coverage_movie( params, map_struct, coverage_struct, plot_sun_moon: bool = True ): - idx = np.isfinite(coverage_struct["data"][:, 4]) if not idx.size: return diff --git a/gwemopt/plotting/plot_efficiency.py b/gwemopt/plotting/plot_efficiency.py index 614051a..f3e8e07 100644 --- a/gwemopt/plotting/plot_efficiency.py +++ b/gwemopt/plotting/plot_efficiency.py @@ -1,4 +1,3 @@ -import healpy as hp import matplotlib.pyplot as plt import numpy as np diff --git a/gwemopt/plotting/plot_skymap.py b/gwemopt/plotting/plot_skymap.py index 16fc979..4ff4111 100644 --- a/gwemopt/plotting/plot_skymap.py +++ b/gwemopt/plotting/plot_skymap.py @@ -1,10 +1,5 @@ -import healpy as hp -import ligo.skymap.plot +import ligo.skymap.plot # noqa: F401 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 def plot_skymap(params, map_struct, colnames=["PROB", "DISTMEAN", "DISTSTD"]): @@ -19,7 +14,6 @@ def plot_skymap(params, map_struct, colnames=["PROB", "DISTMEAN", "DISTSTD"]): for col in colnames: if col in columns: - args = {"projection": params["projection"]} if args["projection"] == "astro globe": args["center"] = map_struct["center"] diff --git a/gwemopt/plotting/plot_tiles.py b/gwemopt/plotting/plot_tiles.py index 10408f2..cd48afc 100644 --- a/gwemopt/plotting/plot_tiles.py +++ b/gwemopt/plotting/plot_tiles.py @@ -1,6 +1,3 @@ -import copy - -import healpy as hp import matplotlib.pyplot as plt import numpy as np from tqdm import tqdm diff --git a/gwemopt/plotting/style.py b/gwemopt/plotting/style.py index 743c5ae..aded6d8 100644 --- a/gwemopt/plotting/style.py +++ b/gwemopt/plotting/style.py @@ -1,8 +1,6 @@ -import healpy as hp -import ligo.skymap.plot +import ligo.skymap.plot # noqa: F401 import matplotlib import matplotlib.pyplot as plt -import numpy as np from astropy.coordinates import get_body from astropy.time import Time diff --git a/gwemopt/run.py b/gwemopt/run.py index b10350c..43f4e7e 100644 --- a/gwemopt/run.py +++ b/gwemopt/run.py @@ -1,4 +1,3 @@ -import os import sys from pathlib import Path @@ -22,7 +21,6 @@ plot_inclination, plot_skymap, ) -from gwemopt.utils import calculate_observability def run(args=None): diff --git a/gwemopt/scheduler.py b/gwemopt/scheduler.py index e59cfc6..c768144 100644 --- a/gwemopt/scheduler.py +++ b/gwemopt/scheduler.py @@ -1,19 +1,14 @@ import copy -import astropy.coordinates -import astropy.units as u import ephem -import ligo.segments as segments import numpy as np from astropy.time import Time from ortools.linear_solver import pywraplp -from gwemopt.tiles import balance_tiles, optimize_max_tiles, schedule_alternating from gwemopt.utils import angular_distance, solve_milp def get_altaz_tile(ra, dec, observer, obstime): - observer.date = ephem.Date(obstime.iso) fxdbdy = ephem.FixedBody() @@ -601,7 +596,6 @@ def scheduler(params, config_struct, tile_struct): config_struct: the telescope configurations tile_struct: the tiles, contains time allocation information """ - import gwemopt.segments # import gwemopt.segments_astroplan coverage_struct = {} @@ -683,7 +677,7 @@ def scheduler(params, config_struct, tile_struct): # calculate airmass for each tile at the start of its exposure: t = Time(mjd_exposure_start, format="mjd") - alt, az = get_altaz_tile( + alt, _ = get_altaz_tile( tile_struct_hold["ra"], tile_struct_hold["dec"], observer, t ) airmass = 1 / np.cos((90.0 - alt) * np.pi / 180) diff --git a/gwemopt/segments.py b/gwemopt/segments.py index 2230915..5f7b747 100644 --- a/gwemopt/segments.py +++ b/gwemopt/segments.py @@ -1,7 +1,5 @@ import copy -import astropy.coordinates -import astropy.units as u import ephem import ligo.segments as segments import numpy as np @@ -11,8 +9,7 @@ from gwemopt.utils.geometry import angular_distance from gwemopt.utils.misc import get_exposures -from gwemopt.utils.parallel import tqdm_joblib -from gwemopt.utils.sidereal_time import greenwich_sidereal_time +from gwemopt.utils.sidereal_time import hour_angle # conversion between MJD (tt) and DJD (what ephem uses) MJD_TO_DJD = -2400000.5 + 2415020.0 @@ -33,7 +30,7 @@ def get_telescope_segments(params): params["config"][telescope]["tot_obs_time"] = 0.0 continue - nexp, junk = np.array(params["config"][telescope]["exposurelist"]).shape + nexp, _ = np.array(params["config"][telescope]["exposurelist"]).shape params["config"][telescope]["n_windows"] = nexp tot_obs_time = ( np.sum(np.diff(np.array(params["config"][telescope]["exposurelist"]))) @@ -111,7 +108,7 @@ def get_moon_segments(config_struct, segmentlist, moon_radecs, radec): return moonsegmentlist -def get_ha_segments(config_struct, segmentlist, observer, radec): +def get_ha_segments(config_struct, segmentlist, radec): if "ha_constraint" in config_struct: ha_constraint = config_struct["ha_constraint"].split(",") ha_min = float(ha_constraint[0]) @@ -131,11 +128,7 @@ def get_ha_segments(config_struct, segmentlist, observer, radec): for seg in segmentlist: mjds = np.linspace(seg[0], seg[1], 100) tt = Time(mjds, format="mjd", scale="utc") - lst = np.mod( - np.deg2rad(config_struct["longitude"]) + greenwich_sidereal_time(tt), - 2 * np.pi, - ) - ha = (lst - np.deg2rad(radec[0])) * 24 / (2 * np.pi) + ha = hour_angle(tt.jd, tt.gps, config_struct["longitude"], radec[0], 0) idx = np.where((ha >= ha_min) & (ha <= ha_max))[0] if len(idx) >= 2: halist.append(segments.segment(mjds[idx[0]], mjds[idx[-1]])) @@ -148,7 +141,6 @@ def get_segments(params, config_struct): event_mjd = Time(gpstime, format="gps", scale="utc").mjd segmentlist = segments.segmentlist() - n_windows = len(params["Tobs"]) // 2 start_segments = event_mjd + params["Tobs"][::2] end_segments = event_mjd + params["Tobs"][1::2] for start_segment, end_segment in zip(start_segments, end_segments): @@ -202,7 +194,6 @@ def get_segments(params, config_struct): def get_segments_tile(config_struct, radec, segmentlist, moon_radecs, airmass): - # check for empty segmentlist and immediately return if len(segmentlist) == 0: return segments.segmentlistdict() @@ -247,7 +238,7 @@ def get_segments_tile(config_struct, radec, segmentlist, moon_radecs, airmass): # moonsegmentlist = get_skybrightness(\ # config_struct,segmentlist,observer,fxdbdy,radec) - halist = get_ha_segments(config_struct, segmentlist, observer, radec) + halist = get_ha_segments(config_struct, segmentlist, radec) moonsegmentlist = get_moon_segments(config_struct, segmentlist, moon_radecs, radec) @@ -292,19 +283,16 @@ def get_segments_tiles(params, config_struct, tile_struct): moon_radecs = get_moon_radecs(segmentlist, observer) if params["doParallel"]: - with tqdm_joblib( - tqdm(desc="segment generation", total=len(radecs)) - ) as progress_bar: - tilesegmentlists = Parallel( - n_jobs=params["Ncores"], - backend=params["parallelBackend"], - batch_size=int(len(radecs) / params["Ncores"]) + 1, - )( - delayed(get_segments_tile)( - config_struct, radec, segmentlist, moon_radecs, params["airmass"] - ) - for radec in radecs + tilesegmentlists = Parallel( + n_jobs=params["Ncores"], + backend=params["parallelBackend"], + batch_size=int(len(radecs) / params["Ncores"]) + 1, + )( + delayed(get_segments_tile)( + config_struct, radec, segmentlist, moon_radecs, params["airmass"] ) + for radec in radecs + ) for ii, key in enumerate(keys): tile_struct[key]["segmentlist"] = tilesegmentlists[ii] else: diff --git a/gwemopt/tests/test_utils.py b/gwemopt/tests/test_utils.py new file mode 100644 index 0000000..65213cd --- /dev/null +++ b/gwemopt/tests/test_utils.py @@ -0,0 +1,51 @@ +import numpy as np +from astropy.time import Time + +from gwemopt.utils import angular_distance, greenwich_sidereal_time, hour_angle + + +def test_greenwich_sidereal_time(): + """ + Test greenwich_sidereal_time function + """ + + # Test with a known date and time + time = Time("2022-01-01T00:00:00", scale="utc") + + gst = greenwich_sidereal_time(time.jd, time.gps, 0) % (2 * np.pi) + + assert np.isclose( + gst, 1.7563325, atol=1e-6 + ), f"Expected gst: 1.7536325, but got: {gst}" + + +def test_hour_angle(): + """ + Test hour_angle function + """ + + # Test with known values + time = Time("2022-01-01T00:00:00", scale="utc") + longitude = 45.0 + ra = 10.0 + ha = hour_angle(time.jd, time.gps, longitude, ra, 0) + + assert np.isclose( + ha, 9.042003, atol=1e-6 + ), f"Expected hour angle: 1.7563325, but got: {ha}" + + +def test_angular_distance(): + """ + Test angular_distance function + """ + + # Test with known coordinates + ra1, dec1 = 0.0, 0.0 # in degrees + ra2, dec2 = 90.0, 0.0 # in degrees + + distance = angular_distance(ra1, dec1, ra2, dec2) + + assert np.isclose( + distance, 90.0, atol=1e-6 + ), f"Expected distance: 90.0, but got: {distance}" diff --git a/gwemopt/tiles.py b/gwemopt/tiles.py index bb228d3..834fe3e 100644 --- a/gwemopt/tiles.py +++ b/gwemopt/tiles.py @@ -1,18 +1,14 @@ -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 ICRS, SkyCoord +from astropy.coordinates import SkyCoord from astropy.time import Time -from joblib import Parallel, delayed from mocpy import MOC from regions import CircleSkyRegion, PolygonSkyRegion, RectangleSkyRegion -from scipy.stats import norm from shapely.geometry import MultiPoint from tqdm import tqdm @@ -20,7 +16,6 @@ import gwemopt.moc import gwemopt.segments from gwemopt.utils.geometry import angular_distance -from gwemopt.utils.parallel import tqdm_joblib TILE_TYPES = ["moc", "galaxy"] @@ -35,10 +30,7 @@ def slice_map_tiles(params, map_struct, coverage_struct): ipix_keep = np.where(csm <= params["iterativeOverlap"])[0] for ii in range(len(coverage_struct["ipix"])): - data = coverage_struct["data"][ii, :] ipix = coverage_struct["ipix"][ii] - FOV = coverage_struct["FOV"][ii] - ipix_slice = np.setdiff1d(ipix, ipix_keep) if len(ipix_slice) == 0: continue @@ -47,7 +39,7 @@ def slice_map_tiles(params, map_struct, coverage_struct): return map_struct -def slice_number_tiles(params, telescope, tile_struct, coverage_struct): +def slice_number_tiles(params, tile_struct, coverage_struct): max_nb_tile = params["max_nb_tiles"] if max_nb_tile is None: return tile_struct, False @@ -82,8 +74,8 @@ def eject_tiles(params, telescope, tile_struct): return tile_struct -def balance_tiles(params, tile_struct, coverage_struct): - filters, exposuretimes = params["filters"], params["exposuretimes"] +def balance_tiles(params, coverage_struct): + filters = params["filters"] keys_scheduled = coverage_struct["data"][:, 5] @@ -194,7 +186,7 @@ def optimize_max_tiles( tile_struct_hold = copy.copy(opt_tile_struct) keys_scheduled = opt_coverage_struct["data"][:, 5] - unique, freq = np.unique(keys_scheduled, return_counts=True) + _, freq = np.unique(keys_scheduled, return_counts=True) n_equal = np.sum(freq == len(params["filters"])) n_dif = np.sum(freq != len(params["filters"])) @@ -232,7 +224,7 @@ def optimize_max_tiles( ) keys_scheduled = coverage_struct_hold["data"][:, 5] - unique, freq = np.unique(keys_scheduled, return_counts=True) + _, freq = np.unique(keys_scheduled, return_counts=True) counter = np.sum(freq == len(params["filters"])) countervals.append(counter) @@ -291,7 +283,7 @@ def optimize_max_tiles( ) keys_scheduled = coverage_struct_hold["data"][:, 5] - unique, freq = np.unique(keys_scheduled, return_counts=True) + _, freq = np.unique(keys_scheduled, return_counts=True) counter = np.sum(freq == len(params["filters"])) countervals.append(counter) @@ -327,9 +319,7 @@ def optimize_max_tiles( tile_struct_hold[key]["prob"] = prob[key] if "epochs" in tile_struct_hold[key]: tile_struct_hold[key]["epochs"] = np.empty((0, 8)) - doReschedule, balanced_fields = balance_tiles( - params_hold, opt_tile_struct, opt_coverage_struct - ) + doReschedule, _ = balance_tiles(params_hold, opt_coverage_struct) params_hold["unbalanced_tiles"] = ( unbalanced_tiles + params_hold["unbalanced_tiles"] ) @@ -347,7 +337,7 @@ def optimize_max_tiles( tile_struct_hold, ) keys_scheduled = coverage_struct["data"][:, 5] - unique, freq = np.unique(keys_scheduled, return_counts=True) + _, freq = np.unique(keys_scheduled, return_counts=True) n_equal = np.sum(freq == len(params["filters"])) n_dif = np.sum(freq != len(params["filters"])) n_1 = n_equal @@ -395,7 +385,7 @@ def check_overlapping_tiles(params, tile_struct, coverage_struct): ra=coverage_ras * u.degree, dec=coverage_decs * u.degree, frame="icrs" ) if params["tilesType"] == "galaxy": - for ii, key in enumerate(keys): + for key in keys: catalog2 = SkyCoord( ra=tile_struct[key]["ra"] * u.degree, dec=tile_struct[key]["dec"] * u.degree, @@ -417,7 +407,7 @@ def check_overlapping_tiles(params, tile_struct, coverage_struct): axis=0, ) else: - for ii, key in enumerate(keys): + for key in keys: catalog2 = SkyCoord( ra=tile_struct[key]["ra"] * u.degree, dec=tile_struct[key]["dec"] * u.degree, @@ -501,7 +491,7 @@ def schedule_alternating( else: filt_change_time = 0 if params["treasuremap_token"] is not None and previous_coverage_struct: - tile_struct_hold = check_overlapping_tiles( + check_overlapping_tiles( params, tile_struct, previous_coverage_struct ) # maps field ids to tile_struct @@ -558,7 +548,7 @@ def schedule_alternating( ) if params["max_nb_tiles"] is not None: tile_struct, doReschedule = slice_number_tiles( - params, telescope, tile_struct, coverage_struct + params, tile_struct, coverage_struct ) if doReschedule: @@ -981,9 +971,7 @@ def compute_tiles_map( func=None, catalog_struct=None, ): - if func is None: - f = lambda x: np.sum(x) - elif func == "center": + if func == "center": keys = tile_struct.keys() ntiles = len(keys) vals = np.nan * np.ones((ntiles,)) @@ -1004,10 +992,6 @@ def compute_tiles_map( val = np.sum(catalog_struct[params["galaxy_grade"]][galaxies]) vals[ii] = val return vals - else: - f = lambda x: eval(func) - - prob = copy.deepcopy(skymap) keys = tile_struct.keys() ntiles = len(keys) diff --git a/gwemopt/utils/__init__.py b/gwemopt/utils/__init__.py index e7aba99..b3357d7 100644 --- a/gwemopt/utils/__init__.py +++ b/gwemopt/utils/__init__.py @@ -4,5 +4,5 @@ from gwemopt.utils.observability import calculate_observability from gwemopt.utils.param_utils import readParamsFromFile from gwemopt.utils.pixels import get_region_moc -from gwemopt.utils.sidereal_time import greenwich_sidereal_time +from gwemopt.utils.sidereal_time import greenwich_sidereal_time, hour_angle from gwemopt.utils.treasuremap import get_treasuremap_pointings diff --git a/gwemopt/utils/geometry.py b/gwemopt/utils/geometry.py index a77dea4..9e3ee93 100644 --- a/gwemopt/utils/geometry.py +++ b/gwemopt/utils/geometry.py @@ -1,6 +1,8 @@ import numpy as np +from numba import njit +@njit def angular_distance(ra1, dec1, ra2, dec2): """ Calculate the angular distance between two points on the sky diff --git a/gwemopt/utils/milp.py b/gwemopt/utils/milp.py index 7ee0022..3c448b7 100644 --- a/gwemopt/utils/milp.py +++ b/gwemopt/utils/milp.py @@ -1,5 +1,3 @@ -import time - import numpy as np import pulp from tqdm import tqdm @@ -16,7 +14,6 @@ def solve_milp( milpSolver="PULP_CBC_CMD", milpOptions={}, ): - cost_matrix_mask = cost_matrix > 10 ** (-10) optimal_points = [] diff --git a/gwemopt/utils/misc.py b/gwemopt/utils/misc.py index 2d1426f..0e48944 100644 --- a/gwemopt/utils/misc.py +++ b/gwemopt/utils/misc.py @@ -1,13 +1,5 @@ -import copy - -import astropy -import astropy.coordinates -import astropy.units as u -import healpy as hp import ligo.segments as segments -import ligo.skymap.distance as ligodist import numpy as np -from astropy.time import Time, TimeDelta def integrationTime(T_obs, pValTiles, func=None, T_int=60.0): diff --git a/gwemopt/utils/parallel.py b/gwemopt/utils/parallel.py deleted file mode 100644 index b51e2f7..0000000 --- a/gwemopt/utils/parallel.py +++ /dev/null @@ -1,22 +0,0 @@ -import contextlib - -import joblib -from tqdm import tqdm - - -@contextlib.contextmanager -def tqdm_joblib(tqdm_object): - """Context manager to patch joblib to report into tqdm progress bar given as argument""" - - class TqdmBatchCompletionCallback(joblib.parallel.BatchCompletionCallBack): - def __call__(self, *args, **kwargs): - tqdm_object.update(n=self.batch_size) - return super().__call__(*args, **kwargs) - - old_batch_callback = joblib.parallel.BatchCompletionCallBack - joblib.parallel.BatchCompletionCallBack = TqdmBatchCompletionCallback - try: - yield tqdm_object - finally: - joblib.parallel.BatchCompletionCallBack = old_batch_callback - tqdm_object.close() diff --git a/gwemopt/utils/param_utils.py b/gwemopt/utils/param_utils.py index 1e39f0a..44d2e1d 100644 --- a/gwemopt/utils/param_utils.py +++ b/gwemopt/utils/param_utils.py @@ -1,11 +1,5 @@ -import glob import os -import astroplan -import astropy -import numpy as np -from astropy import table - def readParamsFromFile(file): """@read gwemopt params file diff --git a/gwemopt/utils/pixels.py b/gwemopt/utils/pixels.py index f4e9557..6529821 100644 --- a/gwemopt/utils/pixels.py +++ b/gwemopt/utils/pixels.py @@ -6,7 +6,6 @@ def get_region_moc(ra, dec, regions, max_depth=12, n_threads=None): - for reg in regions: ra_tmp = reg.vertices.ra dec_tmp = reg.vertices.dec @@ -25,11 +24,6 @@ def get_region_moc(ra, dec, regions, max_depth=12, n_threads=None): stacked = np.stack((ccd_coords.ra.deg, ccd_coords.dec.deg), axis=1) result = stacked.reshape(-1, ccd_coords.ra.deg.shape[1]) lon_lat_list = [row for row in result] - moc = sum( - MOC.from_polygons( - lon_lat_list, max_depth=np.uint8(max_depth), n_threads=n_threads - ) - ) - mocs.append(moc) + mocs.append(sum(MOC.from_polygons(lon_lat_list, False, 10, n_threads))) return mocs diff --git a/gwemopt/utils/sidereal_time.py b/gwemopt/utils/sidereal_time.py index f1289c3..7dcb816 100644 --- a/gwemopt/utils/sidereal_time.py +++ b/gwemopt/utils/sidereal_time.py @@ -1,35 +1,15 @@ -from datetime import datetime, timedelta +from datetime import datetime from math import pi import numpy as np +from numba import njit GPS_EPOCH = datetime(1980, 1, 6, 0, 0, 0) EPOCH_J2000_0_JD = 2451545.0 -_LEAP_SECONDS = np.asarray( - [ - 46828800, - 78364801, - 109900802, - 173059203, - 252028804, - 315187205, - 346723206, - 393984007, - 425520008, - 457056009, - 504489610, - 551750411, - 599184012, - 820108813, - 914803214, - 1025136015, - 1119744016, - 1167264017, - ] -) - - -def greenwich_sidereal_time(tt, equation_of_equinoxes=0): + + +@njit +def greenwich_sidereal_time(jd, gps, equation_of_equinoxes=0): """ Compute the Greenwich mean sidereal time from the GPS time and equation of equinoxes. @@ -38,11 +18,15 @@ def greenwich_sidereal_time(tt, equation_of_equinoxes=0): Parameters ---------- - tt: astropy.time.Time - The astropy time to convert + jd : float + Julian date. + gps : float + GPS time (in seconds). + equation_of_equinoxes : float, optional + Equation of equinoxes (default is 0). """ - t_hi = (tt.jd - EPOCH_J2000_0_JD) / 36525.0 - t_lo = (tt.gps % 1) / (36525.0 * 86400.0) + t_hi = (jd - EPOCH_J2000_0_JD) / 36525.0 + t_lo = (gps % 1) / (36525.0 * 86400.0) t = t_hi + t_lo @@ -55,3 +39,29 @@ def greenwich_sidereal_time(tt, equation_of_equinoxes=0): sidereal_time += 3155760000.0 * t_hi return sidereal_time * pi / 43200.0 + + +@njit +def hour_angle(jd, gps, long, ra, equation_of_equinoxes=0): + """ + Compute the hour angle from the Julian date, GPS time, observer's longitude, + and equation of equinoxes. + + Parameters + ---------- + jd : float + Julian date. + gps : float + GPS time (in seconds). + long : float + Observer's longitude (in degrees). + ra : float + Target's right ascension (in degrees). + equation_of_equinoxes : float, optional + Equation of equinoxes (default is 0). + """ + lst = np.mod( + np.deg2rad(long) + greenwich_sidereal_time(jd, gps, equation_of_equinoxes), + 2 * np.pi, + ) + return (lst - np.deg2rad(ra)) * 24 / (2 * np.pi) diff --git a/pyproject.toml b/pyproject.toml index 78c4ed5..08c7fce 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,6 +29,7 @@ classifiers = [ 'License :: OSI Approved :: GNU General Public License v3 (GPLv3)', ] dependencies = [ + 'numba', 'numpy>=1.7.1', 'scipy>=1.2.1', 'matplotlib>=1.2.0', @@ -54,8 +55,19 @@ dependencies = [ 'tables', 'regions' ] + dynamic = ["version"] +[project.optional-dependencies] +test = [ + 'pytest', + 'pytest-runner', + 'pytest-rerunfailures', + 'coverage', + 'coveralls', + 'freezegun', +] + [project.scripts] gwemopt-run = "gwemopt.run:run"