diff --git a/rubin_sim/maf/maf_contrib/__init__.py b/rubin_sim/maf/maf_contrib/__init__.py index 2c0b18396..2c9235f63 100644 --- a/rubin_sim/maf/maf_contrib/__init__.py +++ b/rubin_sim/maf/maf_contrib/__init__.py @@ -19,7 +19,6 @@ from .tdes_pop_metric import * from .triplet_metric import * from .var_depth_metric import * -from .var_metrics import * from .xrb_metrics import * from .young_stellar_objects_metric import * from .calculate_lsst_field_visibility_astropy import * diff --git a/rubin_sim/maf/maf_contrib/lss_metrics.py b/rubin_sim/maf/maf_contrib/lss_metrics.py index 8504b7c3c..287b1e83b 100644 --- a/rubin_sim/maf/maf_contrib/lss_metrics.py +++ b/rubin_sim/maf/maf_contrib/lss_metrics.py @@ -8,7 +8,8 @@ class GalaxyCountsMetric(BaseMetric): - """Estimate the number of galaxies expected at a particular coadded depth. + """Estimate the number of galaxies expected at a particular (extragalactic) + coadded depth. """ def __init__(self, m5_col="fiveSigmaDepth", nside=128, metric_name="GalaxyCounts", **kwargs): diff --git a/rubin_sim/maf/maf_contrib/lss_obs_strategy/__init__.py b/rubin_sim/maf/maf_contrib/lss_obs_strategy/__init__.py index d8a5dcc0a..f47aa09ba 100644 --- a/rubin_sim/maf/maf_contrib/lss_obs_strategy/__init__.py +++ b/rubin_sim/maf/maf_contrib/lss_obs_strategy/__init__.py @@ -1,8 +1,3 @@ -from .alm_plots import * -from .artificial_structure_calculation import * -from .coadd_m5_analysis import * from .constants_for_pipeline import * from .galaxy_counts_metric_extended import * from .galaxy_counts_with_pixel_calibration import * -from .masking_algorithm_generalized import * -from .os_bias_analysis import * diff --git a/rubin_sim/maf/maf_contrib/lv_dwarfs/lv_dwarfs_metrics.py b/rubin_sim/maf/maf_contrib/lv_dwarfs/lv_dwarfs_metrics.py index e1db87462..eac8c63cf 100644 --- a/rubin_sim/maf/maf_contrib/lv_dwarfs/lv_dwarfs_metrics.py +++ b/rubin_sim/maf/maf_contrib/lv_dwarfs/lv_dwarfs_metrics.py @@ -8,7 +8,6 @@ import os import astropy.units as u -import healpy as hp import numpy as np from astropy.coordinates import SkyCoord from astropy.io import ascii, fits @@ -20,13 +19,14 @@ def generate_known_lv_dwarf_slicer(): - """Read the Karachentsev+ catalog of nearby galaxies, and put the info about them - into a UserPointSlicer object. + """Read the Karachentsev+ catalog of nearby galaxies, + and put the info about them into a UserPointSlicer object. """ filename = os.path.join(get_data_dir(), "maf/lvdwarfs", "lsst_galaxies_1p25to9Mpc_table.fits") lv_dat0 = fits.getdata(filename) - # Keep only galaxies at dec < 35 deg., and with stellar masses > 10^7 M_Sun (and <1e14). + # Keep only galaxies at dec < 35 deg., + # and with stellar masses > 10^7 M_Sun (and <1e14). lv_dat_cuts = (lv_dat0["dec"] < 35.0) & (lv_dat0["MStars"] > 1e7) & (lv_dat0["MStars"] < 1e14) lv_dat = lv_dat0[lv_dat_cuts] @@ -38,11 +38,13 @@ def generate_known_lv_dwarf_slicer(): return slicer -# make a simulated LF for old galaxy of given integrated B, distance modulus mu, in any of filters ugrizY +# make a simulated LF for old galaxy of given integrated B, +# distance modulus mu, in any of filters ugrizY def make__fake_old_galaxy_lf(int_b, mu, filtername): """ - Make a simulated luminosity function for an old (10 Gyr) dwarf galaxy of given - integrated B magnitude, at a given distance modulus, in any of the filters ugrizY. + Make a simulated luminosity function for an old (10 Gyr) dwarf + galaxy of given integrated B magnitude, at a given distance modulus, + in any of the filters ugrizy. Parameters ---------- @@ -56,7 +58,8 @@ def make__fake_old_galaxy_lf(int_b, mu, filtername): if filtername == "y": filtername == "Y" model_bmag = 6.856379 # integrated B mag of the model LF being read - # Read a simulated luminosity function of [M/H]=-1.5, 10 Gyr stellar population: + # Read a simulated luminosity function of [M/H]=-1.5, 10 Gyr stellar + # population: filename = os.path.join(get_data_dir(), "maf/lvdwarfs", "LF_-1.5_10Gyr.dat") LF = ascii.read(filename, header_start=12) mags = LF["magbinc"] @@ -73,8 +76,8 @@ def make__fake_old_galaxy_lf(int_b, mu, filtername): def make_dwarf_lf_dicts(): """ Create dicts containing g- and i-band LFs for simulated dwarfs between - -10 < M_B < +3, so they can simply be looked up rather than having to - recreate them each time. Dict is keyed on M_B value. + -10 < M_B < +3, so they can simply be looked up rather than having to + recreate them each time. Dict is keyed on M_B value. """ lf_dict_i = {} lf_dict_g = {} @@ -96,14 +99,15 @@ def _sum_luminosity(l_fmags, l_fcounts): """ Sum the luminosities from a given luminosity function. - Uses the first bin's magnitude as a reference, sums luminosities relative to - that reference star, then converts back to magnitude at the end. + Uses the first bin's magnitude as a reference, sums luminosities + relative to that reference star, + then converts back to magnitude at the end. Parameters ---------- - l_fmags : np.array, `float` + l_fmags : `np.array,` `float` Magnitude bin values from the simulated LF. - l_fcounts : np.array, `int` + l_fcounts : `np.array`, `int` Number of stars in each magnitude bin. """ magref = l_fmags[0] @@ -156,20 +160,21 @@ def _dwarf_sblimit(glim, ilim, nstars, lf_dict_g, lf_dict_i, distlim, rng): g_l_fmags0, g_l_fcounts0 = lf_dict_g[mbkey] i_l_fcounts = rng.poisson(i_l_fcounts0) g_l_fcounts = rng.poisson(g_l_fcounts0) - i_l_fmags = i_l_fmags0 + distmod_lim # Add the distance modulus to make it apparent mags - g_l_fmags = g_l_fmags0 + distmod_lim # Add the distance modulus to make it apparent mags - # print(i_l_fcounts0-i_l_fcounts) + # Add the distance modulus to make it apparent mags + i_l_fmags = i_l_fmags0 + distmod_lim + # Add the distance modulus to make it apparent mags + g_l_fmags = g_l_fmags0 + distmod_lim gsel = g_l_fmags <= glim isel = i_l_fmags <= ilim ng = np.sum(g_l_fcounts[gsel]) ni = np.sum(i_l_fcounts[isel]) - # print('fake_mb: ',fake_mb, ' ng: ',ng, ' ni: ', ni, ' nstars: ', nstars) fake_mb += 0.1 if fake_mb > -9.9 and (ng > 0) and (ni > 0): gmag_tot = _sum_luminosity(g_l_fmags[gsel], g_l_fcounts[gsel]) - distmod_lim imag_tot = _sum_luminosity(i_l_fmags[isel], i_l_fcounts[isel]) - distmod_lim - # S = m + 2.5logA, where in this case things are in sq. arcmin, so A = 1 arcmin^2 = 3600 arcsec^2 + # S = m + 2.5logA, where in this case things are in sq. arcmin, + # so A = 1 arcmin^2 = 3600 arcsec^2 sbtot_g = distmod_lim + gmag_tot + 2.5 * np.log10(3600.0) sbtot_i = distmod_lim + imag_tot + 2.5 * np.log10(3600.0) mg_lim = gmag_tot @@ -198,35 +203,39 @@ def _dwarf_sblimit(glim, ilim, nstars, lf_dict_g, lf_dict_i, distlim, rng): class LVDwarfsMetric(BaseMetric): """ - Estimate the detection limit in total dwarf luminosity for resolved dwarf galaxies - at a given distance. + Estimate the detection limit in total dwarf luminosity for + resolved dwarf galaxies at a given distance. - This metric class uses simulated luminosity functions of dwarf galaxies with - known (assumed) luminosities to estimate the detection limit (in total dwarf - luminosity, M_V) for resolved dwarf galaxies at a given distance. It can be - applied to either known galaxies with their discrete positions and distances, - or an entire survey simulation with a fixed distance limit. + This metric class uses simulated luminosity functions of dwarf galaxies + with known (assumed) luminosities to estimate the detection limit + (in total dwarf luminosity, M_V) for resolved dwarf galaxies at a + given distance. It can be applied to either known galaxies with their + discrete positions and distances, or an entire survey simulation with + a fixed distance limit. - In the default use (with the KnownLvDwarfsSlicer), it returns detection limits for - a catalog of known local volume dwarfs, from the Karachentsev+ catalog of nearby galaxies. + In the default use (with the KnownLvDwarfsSlicer), + it returns detection limits for a catalog of known local volume dwarfs, + from the Karachentsev+ catalog of nearby galaxies. Parameters ---------- - radius : `float`, default=2.45, - Radius of the field being considered (for discrete fields only). By default, - UserPointSlicer uses a 2.45-deg field radius. - distlim : `float`, - Distance threshold in Mpc for which to calculate the limiting dwarf detection - luminosity. Only needed for healpix slicers, but *required* if healpix is used. - cmd_frac : `float`, default=0.1, - Fraction of the total area of the color-magnitude diagram that is spanned - by the tracer selection criteria. (e.g., the size of a box in color and - magnitude to select RGB-star candidates) - stargal_contamination : `float`, default=0.4, - Fraction of objects in CMD selection region that are actually unresolved - galaxies that were mis-classified as stars. - nsigma : `float`, default=10.0, - Required detection significance to declare a simulated dwarf "detected." + radius : `float` + Radius of the field being considered (for discrete fields only). + By default, UserPointSlicer uses a 2.45-deg field radius. + distlim : `float` + Distance threshold in Mpc for which to calculate the limiting + dwarf detection luminosity. Only needed for healpix slicers, + but *required* if healpix is used. + cmd_frac : `float` + Fraction of the total area of the color-magnitude diagram + that is spanned by the tracer selection criteria. (e.g., + the size of a box in color and magnitude to select RGB-star candidates) + stargal_contamination : `float` + Fraction of objects in CMD selection region that are actually + unresolved galaxies that were mis-classified as stars. + nsigma : `float` + Required detection significance to declare a simulated + dwarf "detected." """ def __init__( @@ -256,7 +265,8 @@ def __init__( self.distlim = None filename = os.path.join(get_data_dir(), "maf/lvdwarfs", "lsst_galaxies_1p25to9Mpc_table.fits") lv_dat0 = fits.getdata(filename) - # Keep only galaxies at dec < 35 deg., and with stellar masses > 10^7 M_Sun. + # Keep only galaxies at dec < 35 deg., + # and with stellar masses > 10^7 M_Sun. lv_dat_cuts = (lv_dat0["dec"] < 35.0) & (lv_dat0["MStars"] > 1e7) & (lv_dat0["MStars"] < 1e14) lv_dat = lv_dat0[lv_dat_cuts] sc_dat = SkyCoord( @@ -281,11 +291,13 @@ def __init__( self.galaxy_counts_metric.scale = 1 cols = [self.m5_col, self.filter_col] - # GalaxyCountsMetric needs the DustMap, and StarDensityMetric needs StellarDensityMap: + # GalaxyCountsMetric needs the DustMap, + # and StarDensityMetric needs StellarDensityMap: maps = ["DustMap", "StellarDensityMap"] super().__init__(col=cols, metric_name=metric_name, maps=maps, units="M_V limit", **kwargs) - # Set up a random number generator, so that metric results are repeatable + # Set up a random number generator, + # so that metric results are repeatable self.rng = np.random.default_rng(seed) def run(self, data_slice, slice_point=None): @@ -296,16 +308,19 @@ def run(self, data_slice, slice_point=None): if len(np.where(gband)[0]) == 0 or len(np.where(iband)[0]) == 0: return self.badval - # calculate the dust-extincted coadded 5-sigma limiting mags in the g and i bands: + # calculate the dust-extincted coadded 5-sigma limiting mags + # in the g and i bands: g5 = self.exgal_coaddm5.run(data_slice[gband], slice_point) i5 = self.exgal_coaddm5.run(data_slice[iband], slice_point) if g5 < 15 or i5 < 15: - # If the limiting magnitudes won't even match the stellar density maps, exit + # If the limiting magnitudes won't even match the + # stellar density maps, exit return self.badval # Find the number of stars per sq arcsecond at the i band limit - # (this is a bit of a hack to get the starDensityMetric to calculate the nstars at this mag exactly) + # (this is a bit of a hack to get the starDensityMetric to + # calculate the nstars at this mag exactly) star_i5 = min(27.9, i5) self.star_density_metric.magLimit = star_i5 @@ -315,21 +330,25 @@ def run(self, data_slice, slice_point=None): ngal_sqdeg = self.galaxy_counts_metric.run(data_slice, slice_point) # GalaxyCountsMetric is undefined in some places. These cases return # zero; catch these and set the galaxy counts in those regions to a - # very high value. (this may not be true after catching earlier no-visits issues) + # very high value. (this may not be true after catching earlier + # no-visits issues) if ngal_sqdeg < 10.0: ngal_sqdeg = 1e7 - # Convert from per sq deg and per sq arcsecond into #'s per sq arcminute + # Convert from per sq deg and per sq arcsecond + # into #'s per sq arcminute ngal_sqarcmin = ngal_sqdeg / 3600 nstar_sqarcmin = nstar_sqarcsec * 3600 if ngal_sqarcmin < 0 or nstar_sqarcmin < 0: print( f"Here be a problem - ngals_sqarcmin {ngal_sqarcmin} or nstar_sqarcmin {nstar_sqarcmin} " - f'are negative. depths: {g5}, {i5}. {slice_point["ra"], slice_point["dec"], slice_point["sid"]}' + f'are negative. depths: {g5}, {i5}. ' + f'{slice_point["ra"], slice_point["dec"], slice_point["sid"]}' ) - # The number of stars required to reach nsigma is nsigma times the Poisson - # fluctuations of the background (stars+galaxies contamination): + # The number of stars required to reach nsigma is + # nsigma times the Poisson fluctuations of the background + # (stars+galaxies contamination): nstars_required = self.nsigma * np.sqrt( (ngal_sqarcmin * self.cmd_frac * self.stargal_contamination) + (nstar_sqarcmin * self.cmd_frac) ) @@ -338,13 +357,12 @@ def run(self, data_slice, slice_point=None): # Use the distlim if a healpix slicer is input distlim = self.distlim else: - # Use discrete distances for known galaxies if a UserPointSlicer: + # Use discrete distances for known galaxies if a + # UserPointSlicer: distlim = slice_point["distance"] * u.Mpc - # sc_slice = SkyCoord(ra=slice_point['ra']*u.rad, dec=slice_point['dec']*u.rad) - # seps = sc_slice.separation(self.sc_dat) - # distlim = self.sc_dat[seps.argmin()].distance - # Calculate the limiting luminosity and surface brightness based on g5 and i5: + # Calculate the limiting luminosity and surface brightness + # based on g5 and i5: mg_lim, mi_lim, sb_g_lim, sb_i_lim, flag_lim = _dwarf_sblimit( g5, i5, @@ -359,8 +377,8 @@ def run(self, data_slice, slice_point=None): mv = self.badval else: - # To go from HSC g and i bands to V, use the conversion from Appendix A - # of Komiyama+2018, ApJ, 853, 29: + # To go from HSC g and i bands to V, use the conversion + # from Appendix A of Komiyama+2018, ApJ, 853, 29: # V = g_hsc - 0.371*(gi_hsc)-0.068 mv = mg_lim - 0.371 * (mg_lim - mi_lim) - 0.068 # sbv = sb_g_lim - 0.371 * (sb_g_lim - sb_i_lim) - 0.068 diff --git a/rubin_sim/maf/maf_contrib/periodic_star_modulation_metric.py b/rubin_sim/maf/maf_contrib/periodic_star_modulation_metric.py index 34a88fa6c..b3c18d6ba 100644 --- a/rubin_sim/maf/maf_contrib/periodic_star_modulation_metric.py +++ b/rubin_sim/maf/maf_contrib/periodic_star_modulation_metric.py @@ -12,38 +12,46 @@ from .periodic_star_metric import PeriodicStar """ This metric is based on the PeriodicStar metric - It was modified in a way to reproduce attempts to identify period/ phase modulation (Blazhko effect) - in RR Lyrae stars. - We are not implementing a period/ phase modulation in the light curve, but rather use short baselines - (e.g.: 20 days) of observations to test how well we can recover the period, phase and amplitude. We - do this as such an attempt is also useful for other purposes, i.e. if we want to test whether we - can just recover period, phase and amplitude from short baselines at all, without necessarily having + It was modified in a way to reproduce attempts to identify + phase modulation (Blazhko effect) in RR Lyrae stars. + We are not implementing a period/ phase modulation in the light curve, + but rather use short baselines (e.g.: 20 days) of observations to test + how well we can recover the period, phase and amplitude. + We do this as such an attempt is also useful for other purposes, + i.e. if we want to test whether we can just recover period, phase + and amplitude from short baselines at all, without necessarily having in mind to look for period/ phase modulations. - Like in the PeriodicStar metric, the light curve of an RR Lyrae star, or a periodic star in general, - is approximated as a simple sin wave. Other solutions might make use of light curve templates to - generate light curves. - Two other modifications we introduced for the PeriodicStarModulationMetric are: - In contrast to the PeriodicStar metric, we allow for a random phase offset to mimic observation - starting at random phase. - Also, we vary the periods and amplitudes within +/- 10 % to allow for a more realistic - sample of variable stars. + Like in the PeriodicStar metric, the light curve of an RR Lyrae star, + or a periodic star in general, is approximated as a simple sin wave. + Other solutions might make use of light curve templates + to generate light curves. + Two other modifications we introduced are: + In contrast to the PeriodicStar metric, we allow for a random phase + offset to mimic observation starting at random phase. + Also, we vary the periods and amplitudes within +/- 10 % to allow + for a more realistic sample of variable stars. This metric is based on the cadence note: - N. Hernitschek, K. Stassun, LSST Cadence Note: Cadence impacts on reliable classification - of standard-candle variable stars (2021) https://docushare.lsst.org/docushare/dsweb/Get/Document-37673 + N. Hernitschek, K. Stassun, LSST Cadence Note: + "Cadence impacts on reliable classification of standard-candle + variable stars (2021)" + https://docushare.lsst.org/docushare/dsweb/Get/Document-37673 """ class PeriodicStarModulationMetric(BaseMetric): - """Evaluate how well a periodic source can be fit on a short baseline, using a Monte Carlo simulation. - - At each slice_point, run a Monte Carlo simulation to see how well a periodic source can be fit. - Assumes a simple sin-wave light-curve, and generates Gaussain noise based in the 5-sigma limiting depth - of each observation. - Light curves are evaluated piecewise to test how well we can recover the period, phase and amplitude - from shorter baselines. We allow for a random phase offset to mimic observation starting at random phase. - Also, we vary the periods and amplitudes within +/- 10 % to allow for a more realistic sample of - variable stars. + """Evaluate how well a periodic source can be fit on a short baseline, + using a Monte Carlo simulation. + + At each slice_point, run a Monte Carlo simulation to see how well a + periodic source can be fit. + Assumes a simple sin-wave light-curve, and generates Gaussain noise + based in the 5-sigma limiting depth of each observation. + Light curves are evaluated piecewise to test how well we can recover + the period, phase and amplitude from shorter baselines. + We allow for a random phase offset to mimic observation starting + at random phase. Also, we vary the periods and amplitudes + within +/- 10 % to allow for a more realistic sample of variable stars. Parameters ---------- @@ -56,11 +64,13 @@ class PeriodicStarModulationMetric(BaseMetric): random_phase : `bool`, opt a random phase is assigned (default False) time_interval : `float`, opt - days (default 50); the interval over which we want to evaluate the light curve + days (default 50); + the interval over which we want to evaluate the light curve n_monte : `int`, opt number of noise realizations to make in the Monte Carlo (default 1000) period_tol : `float`, opt - fractional tolerance on the period to demand for a star to be considered well-fit (default 0.05) + fractional tolerance on the period to demand + for a star to be considered well-fit (default 0.05) amp_tol : `float`, opt fractional tolerance on the amplitude to demand (default 0.10) means : `list` of `float`, opt @@ -120,8 +130,8 @@ def __init__( def run(self, data_slice, slice_point=None): # Bail if we don't have enough points - # (need to fit mean magnitudes in each of the available bands - self.means - # and for a period, amplitude, and phase) + # (need to fit mean magnitudes in each of the available bands + # - self.means and for a period, amplitude, and phase) if data_slice.size < self.means.size + 3: return self.badval @@ -154,8 +164,9 @@ def run(self, data_slice, slice_point=None): mags = self.means + slice_point["distMod"] else: mags = self.means - # slightly different periods and amplitudes (+/- 10 %) to mimic true stars - # random phase offsets to mimic observation starting at random phase + # slightly different periods and amplitudes (+/- 10 %) + # to mimic true stars. random phase offsets to mimic + # observation starting at random phase true_period = random.uniform(0.9, 1.1) * self.period true_amplitude = random.uniform(0.9, 1.1) * self.amplitude if np.isnan(self.phase): @@ -178,7 +189,8 @@ def run(self, data_slice, slice_point=None): fit_obj = PeriodicStar(t_subrun["filter"]) with warnings.catch_warnings(): warnings.simplefilter("ignore") - # If it fails to converge, save values that should fail later + # If it fails to converge, + # save values that should fail later try: parm_vals, pcov = curve_fit( fit_obj, @@ -187,11 +199,12 @@ def run(self, data_slice, slice_point=None): p0=true_params, sigma=dmag, ) - except: + except RuntimeError: parm_vals = true_params * 0 + np.inf fits[i, :] = parm_vals - # Throw out any magnitude fits if there are no observations in that filter + # Throw out any magnitude fits if there are no + # observations in that filter ufilters = np.unique(data_slice[self.filter_col]) if ufilters.size < 9: for key in list(self.filter2index.keys()): diff --git a/rubin_sim/maf/maf_contrib/star_count_mass_metric.py b/rubin_sim/maf/maf_contrib/star_count_mass_metric.py index 1073395b9..7c29011b0 100644 --- a/rubin_sim/maf/maf_contrib/star_count_mass_metric.py +++ b/rubin_sim/maf/maf_contrib/star_count_mass_metric.py @@ -27,6 +27,7 @@ # There are stellar luminosity function maps available within MAF # that may supersede these StarCount functions + class StarCountMassMetric(BaseMetric): """Find the number of stars in a given field in the mass range fainter than magnitude 16 and bright enough to have noise less than @@ -51,7 +52,7 @@ class StarCountMassMetric(BaseMetric): Bandpass to consider. """ - def __init__(self, m1=0.9, m2=1.0, band='i', **kwargs): + def __init__(self, m1=0.9, m2=1.0, band="i", **kwargs): self.m1 = m1 self.m2 = m2 self.band = band diff --git a/rubin_sim/maf/maf_contrib/star_count_metric.py b/rubin_sim/maf/maf_contrib/star_count_metric.py index d7f67cece..4aa158980 100644 --- a/rubin_sim/maf/maf_contrib/star_count_metric.py +++ b/rubin_sim/maf/maf_contrib/star_count_metric.py @@ -25,6 +25,7 @@ # There are stellar luminosity function maps available within MAF # that may supersede these StarCount functions + class StarCountMetric(BaseMetric): """Find the number of stars in a given field between d1 and d2 in parsecs. diff --git a/rubin_sim/maf/maf_contrib/star_counts/coords.py b/rubin_sim/maf/maf_contrib/star_counts/coords.py index 9fa114dd0..0a7c040c9 100644 --- a/rubin_sim/maf/maf_contrib/star_counts/coords.py +++ b/rubin_sim/maf/maf_contrib/star_counts/coords.py @@ -9,7 +9,6 @@ # that uses ephem package, for redundancy purposes. # For use with Field Star Count metric -import sys import numpy as np from scipy.optimize import fsolve @@ -24,23 +23,23 @@ def eq_gal(eq_ra, eq_dec): a = np.radians(eq_ra) def equations(p): - b, l, x = p + b, ll, x = p f1a = np.cos(d) * (np.cos(a - rad1)) f2a = np.sin(d) * np.sin(rad2) + np.cos(d) * np.sin(a - rad1) * np.cos(rad2) f3a = np.sin(d) * np.cos(rad2) - np.cos(d) * np.sin(a - rad1) * np.sin(rad2) - f1 = np.cos(b) * np.cos(l - rad3) - f1a - f2 = np.cos(b) * np.sin(l - rad3) - f2a + f1 = np.cos(b) * np.cos(ll - rad3) - f1a + f2 = np.cos(b) * np.sin(ll - rad3) - f2a f3 = np.sin(b) - f3a return (f1, f2, f3) - b, l, x = fsolve(equations, (0, 0, 0)) + b, ll, x = fsolve(equations, (0, 0, 0)) b_deg = np.degrees(b) % 360 # galactic latitude if b_deg >= 270: b_deg = b_deg - 360 if b_deg > 90: b_deg = 180 - b_deg - l = l + np.pi - l_deg = np.degrees(l) % 360 # galactic longitude + ll = ll + np.pi + l_deg = np.degrees(ll) % 360 # galactic longitude return b_deg, l_deg # http://scienceworld.wolfram.com/astronomy/GalacticCoordinates.html @@ -48,7 +47,6 @@ def equations(p): def eq_gal2(eq_ra, eq_dec): d = np.radians(eq_dec) p = np.radians(eq_ra) - AC = np.radians(90.0) - d AB = np.radians(62.8717) CAB = np.radians(192.8585) - p cos_bc = np.sin(d) * np.cos(AB) + np.cos(d) * np.sin(AB) * np.cos(CAB) @@ -63,9 +61,7 @@ def eq_gal2(eq_ra, eq_dec): cos_cbd = -1 CBD = np.arccos(cos_cbd) b_deg = 90.0 - np.degrees(BC) - ad = np.radians(90.0) cad = np.radians(282.8595) - p - coscd = np.cos(cad) * np.cos(d) coscbd = np.cos(cad) * np.cos(d) / np.sin(BC) if coscbd > 1: coscbd = 1 @@ -100,9 +96,3 @@ def gal_cyn(b_deg, l_deg, dist): rho = np.arctan(y / x) return R, rho, Z - -if __name__ == "__main__": - gal_lat, gal_lon = eq_gal2(float(sys.argv[1]), float(sys.argv[2])) - print(gal_lat, gal_lon) - R, rho, z = gal_cyn(gal_lat, gal_lon, float(sys.argv[3])) - print(R, rho, z) diff --git a/rubin_sim/maf/maf_contrib/star_counts/starcount.py b/rubin_sim/maf/maf_contrib/star_counts/starcount.py index 72a5a394c..c7e8bcdbc 100644 --- a/rubin_sim/maf/maf_contrib/star_counts/starcount.py +++ b/rubin_sim/maf/maf_contrib/star_counts/starcount.py @@ -5,15 +5,10 @@ # Last edited 8/15/2015 # Description: Calculates the number of stars in a given direction and # between a given set of distances. For use with Field Star Count metric -import math -import sys import numpy as np -from scipy.optimize import fsolve - from . import coords, stellardensity -# from rubin_sim.coordUtils import AstronomyBase skyarea = 41253.0 distancebins = 51 @@ -27,21 +22,9 @@ def star_vols(d1, d2, area): def starcount(eq_ra, eq_dec, d1, d2): volumes, distances = star_vols(d1, d2, 9.62) - # b_deg, l_deg=coords.eq_gal2(eq_ra, eq_dec) - # b_deg, l_deg=AstrometryBase.equatorialToGalactic(eq_ra, eq_dec) b_deg, l_deg = coords.eq_gal3(eq_ra, eq_dec) positions = [coords.gal_cyn(b_deg, l_deg, x) for x in distances] densities = [stellardensity.stellardensity(x[0], x[2]) for x in positions] totalcount = np.sum(np.asarray(volumes) * np.asarray(densities)) return totalcount - -if __name__ == "__main__": - print( - starcount( - float(sys.argv[1]), - float(sys.argv[2]), - float(sys.argv[3]), - float(sys.argv[4]), - ) - ) diff --git a/rubin_sim/maf/maf_contrib/star_counts/starcount_bymass.py b/rubin_sim/maf/maf_contrib/star_counts/starcount_bymass.py index 598d30cb9..422b1d1e8 100644 --- a/rubin_sim/maf/maf_contrib/star_counts/starcount_bymass.py +++ b/rubin_sim/maf/maf_contrib/star_counts/starcount_bymass.py @@ -8,12 +8,11 @@ # that will be fainter than mag 16, and have sufficiently low noise # in the given band. For use with Field Star Count metric -import sys import numpy as np -from scipy.optimize import fsolve, minimize, newton +from scipy.optimize import newton -from . import abs_mag, coords, spec_type, stellardensity +from . import abs_mag, spec_type from .starcount import starcount xi = 1.0 @@ -91,15 +90,3 @@ def starcount_bymass(eq_ra, eq_dec, m1, m2, band): distances = [dist_calc(x, band) for x in masses[:-1]] starcounts = [y * starcount(eq_ra, eq_dec, x[0], x[1]) for x, y in zip(distances, massfractions)] return sum(starcounts) - - -if __name__ == "__main__": - print( - starcount_bymass( - float(sys.argv[1]), - float(sys.argv[2]), - float(sys.argv[3]), - float(sys.argv[4]), - sys.argv[5], - ) - ) diff --git a/rubin_sim/maf/maf_contrib/star_counts/stellardensity.py b/rubin_sim/maf/maf_contrib/star_counts/stellardensity.py index d76bce89a..bf3cc0db8 100644 --- a/rubin_sim/maf/maf_contrib/star_counts/stellardensity.py +++ b/rubin_sim/maf/maf_contrib/star_counts/stellardensity.py @@ -5,11 +5,8 @@ # Last edited 8/15/2015 # Description: Calculates the stellar density based off of # Juric et al 2008 and Jackson et al 2002. For use with Field Star Count metric -import math -import sys import numpy as np -from scipy.optimize import fsolve zsun = 25.0 rsun = 8000.0 @@ -58,9 +55,5 @@ def stellardensity(R, Z, rho=0): return tot_density -if __name__ == "__main__": - print(stellardensity(float(sys.argv[1]), float(sys.argv[2]))) - - # Juric et al 2008 # Jackson et al 2002 diff --git a/rubin_sim/maf/maf_contrib/var_depth_metric.py b/rubin_sim/maf/maf_contrib/var_depth_metric.py index 68e6dce4a..15d965b37 100644 --- a/rubin_sim/maf/maf_contrib/var_depth_metric.py +++ b/rubin_sim/maf/maf_contrib/var_depth_metric.py @@ -11,7 +11,23 @@ class VarDepth(BaseMetric): - """Calculate the survey depth that a variable star can be reliably identified.""" + """Calculate the survey depth that a variable star can be + reliably identified. + + Parameters + ---------- + completeness : `float`, opt + Fractional desired completeness of recovered variable sample. + contamination : `float`, opt + Fractional allowed incompleteness of recovered nonvariables. + numruns : `int`, opt + Number of simulated realizations of noise. + Most computationally expensive part of metric. + signal : `float`, opt + Sqrt total pulsational power meant to be recovered. + magres : `float`, opt + desired resolution of variability depth result. + """ def __init__( self, @@ -24,15 +40,6 @@ def __init__( magres=0.01, **kwargs, ): - """ - Instantiate metric. - - :m5col: the column name of the individual visit m5 data. - :completeness: fractional desired completeness of recovered variable sample. - :contamination: fractional allowed incompleteness of recovered nonvariables. - :numruns: number of simulated realizations of noise (most computationally espensive part). - :signal: sqrt total pulsational power meant to be recovered. - :magres: desired resolution of variability depth result.""" self.m5col = m5_col self.completeness = completeness self.contamination = contamination @@ -70,19 +77,20 @@ def run(self, data_slice, slice_point=None): # Since we are treating the underlying signal being representable by a # fixed-width gaussian, its variance pdf is a Chi-squared distribution - # with the degrees of freedom = visits. Since variances add, the variance + # with the degrees of freedom=visits. Since variances add, the variance # pdfs convolve. The cumulative distribution function of the sum of two # random deviates is the convolution of one pdf with a cdf. - # We'll consider the cdf of the noise-only variances because it's easier - # to interpolate + # We'll consider the cdf of the noise-only variances + # because it's easier to interpolate noisesorted = np.sort(noiseonlyvar) # linear interpolation interpnoisecdf = UnivariateSpline( noisesorted, np.arange(self.numruns) / float(self.numruns), k=1, s=0 ) - # We need a binned, signal-only variance probability distribution function for numerical convolution + # We need a binned, signal-only variance probability + # distribution function for numerical convolution numsignalsamples = 100 xsig = np.linspace(chi2.ppf(0.001, N), chi2.ppf(0.999, N), numsignalsamples) signalpdf = chi2.pdf(xsig, N) @@ -90,7 +98,8 @@ def run(self, data_slice, slice_point=None): xsig = (self.signal**2.0) * xsig / N pdfstepsize = xsig[1] - xsig[0] # Since everything is going to use this stepsize down the line, - # normalize so the pdf integrates to 1 when summed (no factor of stepsize needed) + # normalize so the pdf integrates to 1 when summed + # (no factor of stepsize needed) signalpdf /= np.sum(signalpdf) # run through the sample magnitudes, calculate distance between cont @@ -112,11 +121,14 @@ def run(self, data_slice, slice_point=None): # Only do calculation if near the solution: if (len(xnoise) > numsignalsamples / 10) and (not solutionfound): noisecdf = interpnoisecdf(xnoise / scalefact) - noisepdf = noisecdf[1:] - noisecdf[:-1] # turn into a noise pdf + # turn into a noise pdf + noisepdf = noisecdf[1:] - noisecdf[:-1] noisepdf /= np.sum(noisepdf) - xnoise = (xnoise[1:] + xnoise[:-1]) / 2.0 # from cdf to pdf conversion + # from cdf to pdf conversion + xnoise = (xnoise[1:] + xnoise[:-1]) / 2.0 - # calculate and plot the convolution = signal+noise variance dist. + # calculate and plot the convolution = + # signal+noise variance dist. convolution = 0 if len(noisepdf) > len(signalpdf): convolution = np.convolve(noisepdf, signalpdf) diff --git a/rubin_sim/maf/maf_contrib/xrb_metrics.py b/rubin_sim/maf/maf_contrib/xrb_metrics.py index 8a4ec6c9a..3db9796e7 100644 --- a/rubin_sim/maf/maf_contrib/xrb_metrics.py +++ b/rubin_sim/maf/maf_contrib/xrb_metrics.py @@ -24,8 +24,10 @@ def __init__(self, seed=42): def lmxb_abs_mags(self, size=1): """Return LMXB absolute magnitudes per LSST filter. - Absolute magnitude relation is taken from Casares 2018 (2018MNRAS.473.5195C) - Colors are taken from M. Johnson+ 2019 (2019MNRAS.484...19J) + Absolute magnitude relation is taken from Casares 2018 + (2018MNRAS.473.5195C) + Colors are taken from M. Johnson+ 2019 + (2019MNRAS.484...19J) Parameters ---------- @@ -42,7 +44,7 @@ def lmxb_abs_mags(self, size=1): # Derive random orbital periods from the sample in Casares 18 Table 4 # Since there are significant outliers from a single Gaussian sample, - # take random choices with replacement and then perturb them fractionally + # take random choices with replacement, then perturb them fractionally catalog__porb = np.array( [ 33.85, @@ -153,7 +155,8 @@ def fred(self, t, amplitude, tau_rise, tau_decay): return amplitude * np.exp(2 * np.sqrt(tau_rise / tau_decay)) * np.exp(-tau_rise / t - t / tau_decay) def lightcurve(self, t, filtername, params): - """Generate an XRB outburst lightcurve for given times and a single filter. + """Generate an XRB outburst lightcurve for given times + and a single filter. Uses a simple fast-rise, exponential decay with parameters taken from Chen, Shrader, & Livio 1997 (ApJ 491, 312). @@ -173,7 +176,8 @@ def lightcurve(self, t, filtername, params): Returns ------- lc : `array` - Magnitudes of the outburst at the specified times in the given filter + Magnitudes of the outburst at the specified times in + the given filter """ # fill lightcurve with nondetections @@ -194,7 +198,8 @@ def lightcurve(self, t, filtername, params): return lc def detectable_duration(self, params, ebv, distance): - """Determine time range an outburst is detectable with perfect sampling. + """Determine time range an outburst is detectable with + perfect sampling. Does not consider visibility constraints. @@ -210,9 +215,11 @@ def detectable_duration(self, params, ebv, distance): Returns ---------- visible_start_time : `float` - first time relative to outburst start that the outburst could be detected + first time relative to outburst start that the outburst + could be detected visible_end_time : `float` - last time relative to outburst start that the outburst could be detected + last time relative to outburst start that the outburst + could be detected """ nmodelt = 10000 @@ -256,6 +263,24 @@ def detectable_duration(self, params, ebv, distance): class XRBPopMetric(BaseMetric): + """Evaluate whether a given XRB would be detectable. + + Includes a variety of detection criteria options, including if the + XRB is possible to detect, if it is detected at least pts_needed times, + or if it is detected pts_early times within t_early days of the start of + the outburst. + + Parameters + ---------- + pts_needed : `int`, opt + Minimum number of detections, for simple `detected` option. + mjd0 : `float`, opt + Start of survey. + output_lc : `bool`, opt + If True, output lightcurve points. + If False, just return metric values. + """ + def __init__( self, metric_name="XRBPopMetric", @@ -264,6 +289,8 @@ def __init__( filter_col="filter", night_col="night", pts_needed=2, + pts_early=2, + t_early=2, mjd0=None, output_lc=False, badval=-666, @@ -275,6 +302,8 @@ def __init__( self.filter_col = filter_col self.night_col = night_col self.pts_needed = pts_needed + self.pts_early = pts_early + self.t_early = t_early # `bool` variable, if True the light curve will be exported self.output_lc = output_lc @@ -296,7 +325,9 @@ def __init__( self.comment = "Number or characterization of XRBs." def _ever_detect(self, where_detected): - """Simple detection criteria: detect at least a certain number of times""" + """Simple detection criteria: detect at least a certain number + of times. + """ # Detected data points return np.size(where_detected) >= self.pts_needed @@ -324,12 +355,14 @@ def _early_detect(self, where_detected, time, early_window_days=7.0, n_early_det return np.sum(time[where_detected] <= early_window_days) >= n_early_detections def _mean_time_between_detections(self, t): - """Calculate the mean time between detections over the visible interval. + """Calculate the mean time between detections over the + visible interval. Parameters ---------- t : `array` - Times of detections, bracketed by the start and end visibility times + Times of detections, bracketed by the start and + end visibility times Return ---------- @@ -340,7 +373,8 @@ def _mean_time_between_detections(self, t): return np.mean(np.sort(np.diff(t))) def _possible_to_detect(self, visible_duration): - """Return True if the outburst is ever bright enough for LSST to detect + """Return True if the outburst is ever bright enough + for LSST to detect. Parameters ---------- @@ -379,7 +413,7 @@ def run(self, data_slice, slice_point=None): result["possible_to_detect"] = self._possible_to_detect(slice_point["visible_duration"]) result["ever_detect"] = self._ever_detect(where_detected) - result["early_detect"] = self._early_detect(where_detected, t) + result["early_detect"] = self._early_detect(where_detected, t, self.t_early, self.pts_early) result["number_of_detections"] = self._number_of_detections(where_detected) if result["number_of_detections"] > 1: @@ -430,7 +464,7 @@ def reduce_mean_time_between_detections(self, metric): def generate_xrb_pop_slicer(t_start=1, t_end=3652, n_events=10000, seed=42): """Generate a population of XRB events, and put the info about them - into a UserPointSlicer object + into a UserPointSlicer object. Parameters ---------- diff --git a/rubin_sim/maf/maf_contrib/young_stellar_objects_metric.py b/rubin_sim/maf/maf_contrib/young_stellar_objects_metric.py index bf05b1369..de6ad6bb9 100644 --- a/rubin_sim/maf/maf_contrib/young_stellar_objects_metric.py +++ b/rubin_sim/maf/maf_contrib/young_stellar_objects_metric.py @@ -51,31 +51,38 @@ def __call__(self, r): class NYoungStarsMetric(BaseMetric): - """Calculate the distance or number of stars with color uncertainty defined by mags/snrs. + """Calculate the distance or number of stars with + color uncertainty defined by mags/snrs. Parameters ---------- - metric_name : str, opt + metric_name : `str`, opt Default 'young_stars'. - m5_col : str, opt - The default column name for m5 information in the input data. Default fiveSigmaDepth. - filter_col : str, opt + m5_col : `str`, opt + The default column name for m5 information in the input data. + Default fiveSigmaDepth. + filter_col : `str`, opt The column name for the filter information. Default filter. - mags : dict - The absolute magnitude of the object in question. Keys of filter name, values in mags. + mags : `dict`, opt + The absolute magnitude of the object in question. + Keys of filter name, values in mags. Default is for a 0.3 solar mass star at age = 100 Myr. - snrs : dict + snrs : `dict`, opt The SNR to demand for each filter. - galb_limit : float, opt - The galactic latitude above which to return zero (degrees). Default 90. - badval : float, opt - The value to return when the metric value cannot be calculated. Default 0. - return_distance : bool, opt - Whether the metric will return the maximum distance that can be reached for each slice_point, + galb_limit : `float`, opt + The galactic latitude above which to return zero (degrees). + Default 90. + badval : `float`, opt + The value to return when the metric value cannot be calculated. + Default 0. + return_distance : `bool`, opt + Whether the metric will return the maximum distance that + can be reached for each slice_point, or the total number of stars down to mags/snrs. - crowding_error: float, opt + crowding_error : `float`, opt Crowding error that gets passed to CrowdingM5Metric. Default 0.25. - use_2D_extinction: Uses the 2D extinction map instead of the 3D one. Default False. + use_2D_extinction : `bool`, opt + Uses the 2D extinction map instead of the 3D one. Default False. """ def __init__( @@ -109,7 +116,7 @@ def __init__( self.return_distance = return_distance units = "kpc" if self.return_distance else "N stars" super().__init__(cols, metric_name=metric_name, maps=maps, units=units, badval=badval, **kwargs) - # Save R_x values for on-the-fly calculation of dust extinction with map + # Save R_x values for on-the-fly calculation of dust extinction self.r_x = DustValues().r_x.copy() # set return type self.m5_col = m5_col @@ -141,8 +148,8 @@ def run(self, data_slice, slice_point=None): sky_area = np.pi * (np.radians(1.75)) ** 2 # if we are outside the galb_limit, return nothing - # Note we could make this a more complicated function that returns an expected density of - # star forming regions + # Note we could make this a more complicated function that + # returns an expected density of star forming regions if np.abs(slice_point["galb"]) > self.galb_limit: return self.badval @@ -179,7 +186,8 @@ def run(self, data_slice, slice_point=None): filtername=filtername, ) distances.append(dist) - # compute the final distance, limited by whichever filter is most shallow + # compute the final distance, limited by whichever filter is + # most shallow final_distance = np.min(distances, axis=-1) / 1e3 # to kpc if self.return_distance: return final_distance diff --git a/rubin_sim/maf/metric_bundles/metric_bundle.py b/rubin_sim/maf/metric_bundles/metric_bundle.py index 75616eff8..fdf588d06 100644 --- a/rubin_sim/maf/metric_bundles/metric_bundle.py +++ b/rubin_sim/maf/metric_bundles/metric_bundle.py @@ -21,25 +21,33 @@ def create_empty_metric_bundle(): Returns ------- - MetricBundle - An empty metric bundle, configured with just the :class:`BaseMetric` and :class:`BaseSlicer`. + MetricBundle : `MetricBundle` + An empty metric bundle, + configured with just the :class:`BaseMetric` and :class:`BaseSlicer`. """ return MetricBundle(metrics.BaseMetric(), slicers.BaseSlicer(), "") class MetricBundle: - """The MetricBundle is defined by a combination of a (single) metric, slicer and - constraint - together these define a unique combination of an opsim benchmark. - An example would be: a CountMetric, a HealpixSlicer, and a constraint 'filter="r"'. + """Define the "thing" you are measuring, with a combination of + * metric (calculated per data_slice) + * slicer (how to create the data_slices) + * constraint (an optional definition of a large subset of data) - After the metric is evaluated over the slice_points of the slicer, the resulting - metric values are saved in the MetricBundle. + Together these define a unique combination of an opsim benchmark. + An example would be: + a CountMetric, a HealpixSlicer, and a constraint of 'filter="r"'. - The MetricBundle also saves the summary metrics to be used to generate summary - statistics over those metric values, as well as the resulting summary statistic values. + After the metric is evaluated at each slice_point created by the + slicer, the resulting metric values are saved in the MetricBundle. - Plotting parameters and display parameters (for showMaf) are saved in the MetricBundle, - as well as additional info_label such as the opsim run name, and relevant stackers and maps + The MetricBundle also saves the summary metrics to be used + to generate summary statistics over those metric values, + as well as the resulting summary statistic values. + + Plotting parameters and display parameters (for show_maf) are saved + in the MetricBundle, as well as additional info_label such as the + opsim run name, and relevant stackers and maps to apply when calculating the metric values. Parameters @@ -49,33 +57,46 @@ class MetricBundle: slicer : `~rubin_sim.maf.slicer` The Slicer to apply to the incoming visit data (the observations). constraint : `str` or None, opt - A (sql-style) constraint to apply to the visit data, to apply a broad sub-selection. - stacker_list : `list` of `~rubin_sim.maf.stacker`, opt - A list of pre-configured stackers to use to generate additional columns per visit. - These will be generated automatically if needed, but pre-configured versions will override these. + A (sql-style) constraint to apply to the visit data, to apply a + broad sub-selection. + stacker_list : `list` [`~rubin_sim.maf.stacker`], opt + A list of pre-configured stackers to use to generate additional + columns per visit. + These will be generated automatically if needed, but pre-configured + versions will override these. run_name : `str`, opt - The name of the simulation being run. This will be added to output files and plots. - Setting it prevents file conflicts when running the same metric on multiple simulations, and + The name of the simulation being run. + This will be added to output files and plots. + Setting it prevents file conflicts when running the same + metric on multiple simulations, and provides a way to identify which simulation is being analyzed. metadata : `str`, opt A deprecated version of info_label (below). - Values set by metadata will be used for info_label. If both are set, info_label is used. + Values set by metadata will be used for info_label. + If both are set, info_label is used. info_label : `str` or None, opt Information to add to the output metric data file name and plot labels. - If this is not provided, it will be auto-generated from the constraint (if any). - Setting this provides an easy way to specify different configurations of a metric, a slicer, + If this is not provided, it will be auto-generated from the + constraint (if any). + Setting this provides an easy way to specify different + configurations of a metric, a slicer, or just to rewrite your constraint into friendlier terms. - (i.e. a constraint like 'note not like "%DD%"' can become "non-DD" in the file name and plot labels + (i.e. a constraint like 'note not like "%DD%"' can become + "non-DD" in the file name and plot labels by specifying info_label). plot_dict : `dict` of plotting parameters, opt Specify general plotting parameters, such as x/y/color limits. display_dict : `dict` of display parameters, opt - Specify parameters for showMaf web pages, such as the side bar labels and figure captions. - Keys: 'group', 'subgroup', 'caption', and 'order' (such as to set metrics in filter order, etc) + Specify parameters for show_maf web pages, such as the + side bar labels and figure captions. + Keys: 'group', 'subgroup', 'caption', and 'order' + (such as to set metrics in filter order, etc) summary_metrics : `list` of `~rubin_sim.maf.metrics` - A list of summary metrics to run to summarize the primary metric, such as MedianMetric, etc. + A list of summary metrics to run to summarize the + primary metric, such as MedianMetric, etc. maps_list : `list` of `~rubin_sim.maf.maps` - A list of pre-configured maps to use for the metric. This will be auto-generated if specified + A list of pre-configured maps to use for the metric. + This will be auto-generated if specified by the metric class, but pre-configured versions will override these. """ @@ -144,7 +165,7 @@ def __init__( map_names = [map_name.__class__.__name__ for map_name in self.maps_list] if hasattr(self.metric, "maps"): for map_needed in self.metric.maps: - if type(map_needed) == str: + if isinstance(map_needed, str): if map_needed not in map_names: temp_map = getattr(maps, map_needed)() self.maps_list.append(temp_map) @@ -215,7 +236,8 @@ def _setup_metric_values(self): def _build_metadata(self, info_label, metadata=None): """If no info_label is provided, process the constraint - (by removing extra spaces, quotes, the word 'filter' and equal signs) to make a info_label version. + (by removing extra spaces, quotes, the word 'filter' and equal signs) + to make a info_label version. e.g. 'filter = "r"' becomes 'r' """ # Pass the deprecated version into info_label if info_label is not set @@ -234,7 +256,8 @@ def _build_metadata(self, info_label, metadata=None): def _build_file_root(self): """ - Build an auto-generated output filename root (i.e. minus the plot type or .npz ending). + Build an auto-generated output filename root + (i.e. minus the plot type or .npz ending). """ # Build basic version. self.file_root = "_".join( @@ -250,10 +273,10 @@ def _build_file_root(self): def _find_req_cols(self): """Find the columns needed by the metrics, slicers, and stackers. - If there are any additional stackers required, instatiate them and add them to - the self.stackers list. - (default stackers have to be instantiated to determine what additional columns - are needed from database). + If there are any additional stackers required, instatiate them + and add them to the self.stackers list. + (default stackers have to be instantiated to determine + what additional columns are needed from database). """ # Find all the columns needed by metric and slicer. known_cols = self.slicer.columns_needed + list(self.metric.col_name_arr) @@ -268,7 +291,8 @@ def _find_req_cols(self): if self.col_info.get_data_source(col) == self.col_info.default_data_source: self.db_cols.add(col) else: - # New default stackers could come from metric/slicer or stackers. + # New default stackers could come from metric/slicer + # or stackers. new_stackers.add(self.col_info.get_data_source(col)) # Remove already-specified stackers from default list. for s in self.stacker_list: @@ -304,7 +328,7 @@ def set_summary_metrics(self, summary_metrics): Parameters ---------- - summary_metrics : `List` of [`BaseMetric`] + summary_metrics : `List` [`BaseMetric`] Instantiated summary metrics to use to calculate summary statistics for this metric. """ @@ -448,9 +472,9 @@ def set_run_name(self, run_name, update_file_root=True): Parameters ---------- - run_name: `str` + run_name : `str` Run Name, which will become part of the fileRoot. - fileRoot: `bool`, optional + fileRoot : `bool`, optional Flag to update the fileRoot with the run_name. """ self.run_name = run_name @@ -508,7 +532,8 @@ def write(self, comment="", out_dir=".", outfile_suffix=None, results_db=None): self.write_db(results_db=results_db) def output_json(self): - """Set up and call the baseSlicer outputJSON method, to output to IO string. + """Set up and call the baseSlicer outputJSON method, + to output to IO string. Returns ------- @@ -581,7 +606,7 @@ def load(cls, filename): Parameters ---------- - filename : str + filename : `str` The file from which to read the metric bundle data. """ metric_bundle = cls(metrics.BaseMetric(), slicers.BaseSlicer(), "") @@ -589,28 +614,32 @@ def load(cls, filename): return metric_bundle def compute_summary_stats(self, results_db=None): - """Compute summary statistics on metric_values, using summaryMetrics (metricbundle list). + """Compute summary statistics on metric_values, + using summaryMetrics (metricbundle list). Parameters ---------- results_db : Optional[ResultsDb] - ResultsDb object to use to store the summary statistic values on disk. + ResultsDb object to use to store the + summary statistic values on disk. """ if self.summary_values is None: self.summary_values = {} if self.summary_metrics is not None: - # Build array of metric values, to use for (most) summary statistics. + # Build array of metric values, to use for summary statistics. rarr_std = np.array( list(zip(self.metric_values.compressed())), dtype=[("metricdata", self.metric_values.dtype)], ) for m in self.summary_metrics: - # The summary metric colname should already be set to 'metricdata', but in case it's not: + # The summary metric colname should already be set + # to 'metricdata', but in case it's not: m.colname = "metricdata" summary_name = m.name.replace(" metricdata", "").replace(" None", "") if hasattr(m, "mask_val"): - # summary metric requests to use the mask value, as specified by itself, - # rather than skipping masked vals. + # summary metric requests to use the mask value, + # as specified by itself, + # rather than skipping masked vals. rarr = np.array( list(zip(self.metric_values.filled(m.mask_val))), dtype=[("metricdata", self.metric_values.dtype)], @@ -644,25 +673,28 @@ def reduce_metric( reduce_display_dict=None, ): """Run 'reduceFunc' (any function that operates on self.metric_values). - Typically reduceFunc will be the metric reduce functions, as they are tailored to expect the - metric_values format. - reduceDisplayDict and reducePlotDicts are displayDicts and plotDicts to be - applied to the new metricBundle. + + Typically reduceFunc will be the metric reduce functions, + as they are tailored to expect the metric_values format. + reduceDisplayDict and reducePlotDicts are displayDicts + and plotDicts to be applied to the new metricBundle. Parameters ---------- - reduce_func : Func - Any function that will operate on self.metric_values (typically metric.reduce* function). - reduce_plot_dict : Optional[dict] + reduce_func : `Func` + Any function that will operate on self.metric_values + (typically metric.reduce* function). + reduce_plot_dict : `dict`, opt Plot dictionary for the results of the reduce function. - reduce_display_dict : Optional[dict] + reduce_display_dict : `dict`, opt Display dictionary for the results of the reduce function. Returns ------- - MetricBundle - New metric bundle, inheriting info_label from this metric bundle, but containing the new - metric values calculated with the 'reduceFunc'. + newmetric_bundle: `MetricBundle` + New metric bundle, inheriting info_label from this metric bundle, + but containing the new metric values calculated with + the 'reduceFunc'. """ # Generate a name for the metric values processed by the reduceFunc. if reduce_func_name is not None: @@ -670,7 +702,8 @@ def reduce_metric( else: r_name = reduce_func.__name__.replace("reduce_", "") reduce_name = self.metric.name + "_" + r_name - # Set up metricBundle to store new metric values, and add plot_dict/display_dict. + # Set up metricBundle to store new metric values, + # and add plot_dict/display_dict. newmetric = deepcopy(self.metric) newmetric.name = reduce_name newmetric.metric_dtype = "float" @@ -693,22 +726,27 @@ def reduce_metric( ) # Build a new output file root name. newmetric_bundle._build_file_root() - # Add existing plot_dict (except for title/xlabels etc) into new plot_dict. + # Add existing plot_dict (except for title/xlabels etc) + # into new plot_dict. for k, v in self.plot_dict.items(): if k not in newmetric_bundle.plot_dict: newmetric_bundle.plot_dict[k] = v - # Update newmetric_bundle's plot dictionary with any set explicitly by reducePlotDict. + # Update newmetric_bundle's plot dictionary with + # any set explicitly by reducePlotDict. newmetric_bundle.set_plot_dict(reduce_plot_dict) # Copy the parent metric's display dict into the reduce display dict. newmetric_bundle.set_display_dict(self.display_dict) - # Set the reduce function display 'order' (this is set in the BaseMetric - # by default, but can be overriden in a metric). + # Set the reduce function display 'order' + # (this is set in the BaseMetric by default, + # but can be overriden in a metric). order = newmetric.reduce_order[r_name] newmetric_bundle.display_dict["order"] = order - # And then update the newmetric_bundle's display dictionary with any set + # And then update the newmetric_bundle's + # display dictionary with any set # explicitly by reduceDisplayDict. newmetric_bundle.set_display_dict(reduce_display_dict) - # Set up new metricBundle's metric_values masked arrays, copying metricValue's mask. + # Set up new metricBundle's metric_values masked arrays, + # copying metricValue's mask. newmetric_bundle.metric_values = ma.MaskedArray( data=np.empty(len(self.slicer), "float"), mask=self.metric_values.mask.copy(), @@ -726,24 +764,29 @@ def reduce_metric( def plot(self, plot_handler=None, plot_func=None, outfile_suffix=None, savefig=False): """ - Create all plots available from the slicer. plotHandler holds the output directory info, etc. + Create all plots available from the slicer. + plotHandler holds the output directory info, etc. Parameters ---------- - plot_handler : Optional[PlotHandler] - The plot_handler saves the output location and results_db connection for a set of plots. - plot_func : Optional[BasePlotter] - Any plotter function. If not specified, the plotters in self.plotFuncs will be used. - outfile_suffix : Optional[str] + plot_handler : `~maf.plots.plot_handler`, opt + The plot_handler saves the output location + and results_db connection for a set of plots. + plot_func : `maf.plots.BasePlotter`, opt + Any plotter function. + If not specified, the plotters in self.plotFuncs will be used. + outfile_suffix : `str`, opt Optional string to append to the end of the plot output files. Useful when creating sequences of images for movies. - savefig : Optional[bool] - Flag indicating whether or not to save the figure to disk. Default is False. + savefig : `bool`, opt + Flag indicating whether or not to save the figure to disk. + Default is False. Returns ------- - dict - Dictionary of plot_type:figure number key/value pairs, indicating what plots were created + made_plots : `dict` + Dictionary of plot_type:figure number key/value pairs, + indicating what plots were created and what matplotlib figure numbers were used. """ # Generate a plot_handler if none was set. diff --git a/rubin_sim/maf/metric_bundles/metric_bundle_group.py b/rubin_sim/maf/metric_bundles/metric_bundle_group.py index 40945ecdf..a22923d0b 100644 --- a/rubin_sim/maf/metric_bundles/metric_bundle_group.py +++ b/rubin_sim/maf/metric_bundles/metric_bundle_group.py @@ -18,14 +18,15 @@ def make_bundles_dict_from_list(bundle_list): - """Utility to convert a list of MetricBundles into a dictionary, keyed by the fileRoot names. + """Utility to convert a list of MetricBundles into a dictionary, + keyed by the fileRoot names. Raises an exception if the fileroot duplicates another metricBundle. (Note this should alert to potential cases of filename duplication). Parameters ---------- - bundle_list : `list` of `MetricBundles` + bundle_list : `list` [`MetricBundles`] """ b_dict = {} for b in bundle_list: @@ -36,47 +37,57 @@ def make_bundles_dict_from_list(bundle_list): class MetricBundleGroup: - """The MetricBundleGroup exists to calculate the metric values for a group of - MetricBundles. + """Calculate all values for a group of MetricBundles. - The MetricBundleGroup will query data from a single database table (for multiple - constraints), use that data to calculate metric values for multiple slicers, - and calculate summary statistics and generate plots for all metrics included in + The MetricBundleGroup will query data from a single database table + (for multiple constraints), use that data to calculate metric values + for multiple slicers, and calculate summary statistics and + generate plots for all metrics included in the dictionary passed to the MetricBundleGroup. - We calculate the metric values here, rather than in the individual MetricBundles, - because it is much more efficient to step through a slicer once (and calculate all - the relevant metric values at each point) than it is to repeat this process multiple times. + We calculate the metric values here, rather than in the + individual MetricBundles, because it is much more efficient to step + through a slicer once (and calculate all the relevant metric values + at each point) than it is to repeat this process multiple times. - The MetricBundleGroup also determines how to efficiently group the MetricBundles - to reduce the number of sql queries of the database, grabbing larger chunks of data at once. + The MetricBundleGroup also determines how to efficiently group + the MetricBundles to reduce the number of sql queries of the database, + grabbing larger chunks of data at once. Parameters ---------- - bundle_dict : `dict` or `list` of `MetricBundles` - Individual MetricBundles should be placed into a dictionary, and then passed to - the MetricBundleGroup. The dictionary keys can then be used to identify MetricBundles - if needed -- and to identify new MetricBundles which could be created if 'reduce' - functions are run on a particular MetricBundle. - A bundle_dict can be conveniently created from a list of MetricBundles using - makeBundlesDictFromList (done automatically if a list is passed in) + bundle_dict : `dict` or `list` [`MetricBundles`] + Individual MetricBundles should be placed into a dictionary, + and then passed to the MetricBundleGroup. + The dictionary keys can then be used to identify MetricBundles + if needed -- and to identify new MetricBundles which could be + created if 'reduce' functions are run on a particular MetricBundle. + A bundle_dict can be conveniently created from a list of MetricBundles + using makeBundlesDictFromList (done automatically if a list is passed). db_con : `str` or database connection object - A str that is the path to a sqlite3 file or a database object that can be used by pandas.read_sql. - Advanced use: It is possible to set this to None, in which case data should be passed - directly to the runCurrent method (and runAll should not be used). - out_dir : `str`, optional + A str that is the path to a sqlite3 file or a database object + that can be used by pandas.read_sql. + Advanced use: It is possible to set this to None, in which case + data should be passed directly to the runCurrent method + (and runAll should not be used). + out_dir : `str`, opt Directory to save the metric results. Default is the current directory. - results_db : `ResultsDb`, optional - A results database. If not specified, one will be created in the out_dir. - This database saves information about the metrics calculated, including their summary statistics. - verbose : `bool`, optional + results_db : `ResultsDb`, opt + A results database to store summary stat information. + If not specified, one will be created in the out_dir. + This database saves information about the metrics calculated, + including their summary statistics. + verbose : `bool`, opt Flag to turn on/off verbose feedback. - save_early : `bool`, optional - If True, metric values will be saved immediately after they are first calculated (to prevent - data loss) as well as after summary statistics are calculated. - If False, metric values will only be saved after summary statistics are calculated. - db_table : `str`, optional + save_early : `bool`, opt + If True, metric values will be saved immediately after + they are first calculated (to prevent data loss) as well as after + summary statistics are calculated. + If False, metric values will only be saved after summary statistics + are calculated. + db_table : `str`, opt The name of the table in the db_obj to query for data. + For modern opsim outputs, this table is `observations` (default None). """ def __init__( @@ -90,7 +101,7 @@ def __init__( db_table=None, ): """Set up the MetricBundleGroup.""" - if type(bundle_dict) is list: + if isinstance(bundle_dict, list): bundle_dict = make_bundles_dict_from_list(bundle_dict) # Print occasional messages to screen. self.verbose = verbose @@ -129,19 +140,21 @@ def __init__( def _check_compatible(self, metric_bundle1, metric_bundle2): """Check if two MetricBundles are "compatible". - Compatible indicates that the sql constraints, the slicers, and the maps are the same, and + + Compatible indicates that the sql constraints, the slicers, + and the maps are the same, and that the stackers do not interfere with each other (i.e. are not trying to set the same column in different ways). Returns True if the MetricBundles are compatible, False if not. Parameters ---------- - metric_bundle1 : MetricBundle - metric_bundle2 : MetricBundle + metric_bundle1 : `MetricBundle` + metric_bundle2 : `MetricBundle` Returns ------- - bool + match : `bool` """ if metric_bundle1.constraint != metric_bundle2.constraint: return False @@ -151,7 +164,8 @@ def _check_compatible(self, metric_bundle1, metric_bundle2): return False for stacker in metric_bundle1.stacker_list: for stacker2 in metric_bundle2.stacker_list: - # If the stackers have different names, that's OK, and if they are identical, that's ok. + # If the stackers have different names, that's OK, + # and if they are identical, that's ok. if (stacker.__class__.__name__ == stacker2.__class__.__name__) & (stacker != stacker2): return False # But if we got this far, everything matches. @@ -160,8 +174,8 @@ def _check_compatible(self, metric_bundle1, metric_bundle2): def _find_compatible_lists(self): """Find sets of compatible metricBundles from the currentBundleDict.""" # CompatibleLists stores a list of lists; - # each (nested) list contains the bundleDict _keys_ of a compatible set of metricBundles. - # + # each (nested) list contains the bundleDict _keys_ + # of a compatible set of metricBundles. compatible_lists = [] for k, b in self.current_bundle_dict.items(): found_compatible = False @@ -169,16 +183,20 @@ def _find_compatible_lists(self): comparison_metric_bundle_key = compatible_list[0] compatible = self._check_compatible(self.bundle_dict[comparison_metric_bundle_key], b) if compatible: - # Must compare all metricBundles in each subset (if they are a potential match), - # as the stackers could be different (and one could be incompatible, + # Must compare all metricBundles in each subset + # (if they are a potential match), + # as the stackers could be different + # (and one could be incompatible, # not necessarily the first) for comparison_metric_bundle_key in compatible_list[1:]: compatible = self._check_compatible(self.bundle_dict[comparison_metric_bundle_key], b) if not compatible: - # If we find one which is not compatible, stop and go on to the + # If we find one which is not compatible, + # stop and go on to the # next subset list. break - # Otherwise, we reached the end of the subset and they were all compatible. + # Otherwise, we reached the end of the subset + # and they were all compatible. found_compatible = True compatible_list.append(k) if not found_compatible: @@ -191,23 +209,25 @@ def _find_compatible_lists(self): self.compatible_lists = compatible_lists def run_all(self, clear_memory=False, plot_now=False, plot_kwargs=None): - """Runs all the metricBundles in the metricBundleGroup, over all constraints. + """Runs all the metricBundles in the metricBundleGroup, + over all constraints. - Calculates metric values, then runs reduce functions and summary statistics for - all MetricBundles. + Calculates metric values, then runs reduce functions and summary + statistics for all MetricBundles. Parameters ---------- clear_memory : `bool`, optional - If True, deletes metric values from memory after running each constraint group. + If True, deletes metric values from memory after running + each constraint group. plot_now : `bool`, optional If True, plots the metric values immediately after calculation. plot_kwargs : `bool`, optional kwargs to pass to plotCurrent. """ for constraint in self.constraints: - # Set the 'currentBundleDict' which is a dictionary of the metricBundles which match this - # constraint. + # Set the 'currentBundleDict' which is a dictionary of the + # metricBundles which match this constraint. self.run_current( constraint, clear_memory=clear_memory, @@ -216,13 +236,15 @@ def run_all(self, clear_memory=False, plot_now=False, plot_kwargs=None): ) def set_current(self, constraint): - """Utility to set the currentBundleDict (i.e. a set of metricBundles with the same SQL constraint). + """Utility to set the currentBundleDict + (i.e. a set of metricBundles with the same SQL constraint). Parameters ---------- constraint : `str` - The subset of MetricBundles with metricBundle.constraint == constraint will be - included in a subset identified as the currentBundleDict. + The subset of MetricBundles with metricBundle.constraint == + constraint will be included in a subset identified as the + currentBundleDict. These are the active metrics to be calculated and plotted, etc. """ if constraint is None: @@ -256,17 +278,17 @@ def run_current( ---------- constraint : `str` constraint to use to set the currently active metrics - sim_data : `numpy.ndarray`, optional + sim_data : `np.ndarray`, ops If simData is not None, then this numpy structured array is used instead of querying data from the dbObj. - clear_memory : `bool`, optional + clear_memory : `bool`, ops If True, metric values are deleted from memory after they are calculated (and saved to disk). - plot_now : `bool`, optional + plot_now : `bool`, ops Plot immediately after calculating metric values (instead of the usual procedure, which is to plot after metric values are calculated for all constraints). - plot_kwargs : kwargs, optional + plot_kwargs : kwargs, ops Plotting kwargs to pass to plotCurrent. """ self.set_current(constraint) @@ -371,10 +393,11 @@ def _run_compatible(self, compatible_list): """Runs a set of 'compatible' metricbundles in the MetricBundleGroup dictionary identified by 'compatible_list' keys. - A compatible list of MetricBundles is a subset of the currentBundleDict. + A compatible list of MetricBundles is a subset of the + currentBundleDict. The currentBundleDict == set of MetricBundles with the same constraint. The compatibleBundles == set of MetricBundles with the same constraint, - the same slicer, the same maps applied to the slicer, + AND the same slicer, the same maps applied to the slicer, and stackers which do not clobber each other's data. This is where the work of calculating the metric values is done. @@ -541,12 +564,15 @@ def reduce_current(self, update_summaries=True): # Create a temporary dictionary to hold the reduced metricbundles. reduce_bundle_dict = {} for b in self.current_bundle_dict.values(): - # If there are no reduce functions associated with the metric, skip this metricBundle. + # If there are no reduce functions associated with the metric, + # skip this metricBundle. if len(b.metric.reduce_funcs) > 0: - # Apply reduce functions, creating a new metricBundle in the process (new metric values). + # Apply reduce functions, creating a new metricBundle in + # the process (new metric values). for r in b.metric.reduce_funcs: newmetricbundle = b.reduce_metric(b.metric.reduce_funcs[r], reduce_func_name=r) - # Add the new metricBundle to our metricBundleGroup dictionary. + # Add the new metricBundle to our metricBundleGroup + # dictionary. name = newmetricbundle.metric.name if name in self.bundle_dict: name = newmetricbundle.file_root @@ -560,7 +586,8 @@ def reduce_current(self, update_summaries=True): b.summary_metrics = [] # Add the new metricBundles to the MetricBundleGroup dictionary. self.bundle_dict.update(reduce_bundle_dict) - # And add to to the currentBundleDict too, so we run as part of 'summaryCurrent'. + # And add to to the currentBundleDict too, so we run as part + # of 'summaryCurrent'. self.current_bundle_dict.update(reduce_bundle_dict) def summary_all(self): @@ -575,7 +602,8 @@ def summary_all(self): def summary_current(self): """Run summary statistics on all the metricBundles in the - currently active set of MetricBundles.""" + currently active set of MetricBundles. + """ for b in self.current_bundle_dict.values(): b.compute_summary_stats(self.results_db) @@ -692,7 +720,8 @@ def plot_current( def write_all(self): """Save all the MetricBundles to disk. - Saving all MetricBundles to disk at this point assumes that clearMemory was False. + Saving all MetricBundles to disk at this point assumes that + clearMemory was False. """ for constraint in self.constraints: self.set_current(constraint) diff --git a/rubin_sim/maf/metric_bundles/mo_metric_bundle.py b/rubin_sim/maf/metric_bundles/mo_metric_bundle.py index 3dba0b7c2..17f9c8084 100644 --- a/rubin_sim/maf/metric_bundles/mo_metric_bundle.py +++ b/rubin_sim/maf/metric_bundles/mo_metric_bundle.py @@ -26,36 +26,41 @@ def create_empty_mo_metric_bundle(): Returns ------- ~rubin_sim.maf.metricBundles.MoMetricBundle - An empty metric bundle, configured with just the :class:`BaseMetric` and :class:`BaseSlicer`. + An empty metric bundle, configured with just + the :class:`BaseMetric` and :class:`BaseSlicer`. """ return MoMetricBundle(BaseMoMetric(), MoObjSlicer(), None) def make_completeness_bundle(bundle, completeness_metric, h_mark=None, results_db=None): - """ - Make a mock metric bundle from a bundle which had MoCompleteness or MoCumulativeCompleteness summary - metrics run. This lets us use the plotHandler + plots.MetricVsH to generate plots. - Will also work with completeness metric run in order to calculate fraction of the population, - or with MoCompletenessAtTime metric. + """Make a mock metric bundle from a bundle which had + MoCompleteness or MoCumulativeCompleteness summary + metrics run. This lets us use the plotHandler + plots.MetricVsH + to generate plots. + Will also work with completeness metric run in order to calculate + fraction of the population, or with MoCompletenessAtTime metric. Parameters ---------- - bundle : ~rubin_sim.maf.metricBundles.MetricBundle + bundle : `~rubin_sim.maf.metricBundles.MetricBundle` The metric bundle with a completeness summary statistic. - completeness_metric : ~rubin_sim.maf.metric + completeness_metric : `~rubin_sim.maf.metric` The summary (completeness) metric to run on the bundle. - h_mark : float, optional - The Hmark value to add to the plotting dictionary of the new mock bundle. Default None. - results_db : ~rubin_sim.maf.db.ResultsDb, optional - The results_db in which to record the summary statistic value at Hmark. Default None. + h_mark : `float`, optional + The Hmark value to add to the plotting dictionary of the new + mock bundle. Default None. + results_db : `~rubin_sim.maf.db.ResultsDb`, optional + The results_db in which to record the summary statistic value at + Hmark. Default None. Returns ------- - ~rubin_sim.maf.metricBundles.MoMetricBundle + mo_metric_bundle : `~rubin_sim.maf.metricBundles.MoMetricBundle` """ bundle.set_summary_metrics(completeness_metric) - # This step adds summary values at each point to the original metric - we use this to populate - # the completeness values in the next step. However, we may not want them to go into the results_db. + # This step adds summary values at each point to the original metric - + # we use this to populate the completeness values in the next step. + # However, we may not want them to go into the results_db. bundle.compute_summary_stats(results_db) summary_name = completeness_metric.name # Make up the bundle, including the metric values. @@ -109,9 +114,6 @@ def __init__( child_metrics=None, summary_metrics=None, ): - """ - Instantiate moving object metric bundle, save metric/slicer/constraint, etc. - """ self.metric = metric self.slicer = slicer if constraint == "": @@ -187,7 +189,8 @@ def _reset_metric_bundle(self): self.summary_values = None def _build_metadata(self, info_label): - """If no info_label is provided, auto-generate it from the obs_file + constraint.""" + """If no info_label is provided, auto-generate it from the + obs_file + constraint.""" if info_label is None: try: self.info_label = self.slicer.obsfile.replace(".txt", "").replace(".dat", "") @@ -207,7 +210,8 @@ def _find_req_cols(self): def set_child_bundles(self, child_metrics=None): """ Identify any child metrics to be run on this (parent) bundle. - and create the new metric bundles that will hold the child values, linking to this bundle. + and create the new metric bundles that will hold the child values, + linking to this bundle. Remove the summaryMetrics from self afterwards. """ self.child_bundles = {} @@ -232,12 +236,13 @@ def set_child_bundles(self, child_metrics=None): def compute_summary_stats(self, results_db=None): """ - Compute summary statistics on metric_values, using summaryMetrics, for self and child bundles. + Compute summary statistics on metric_values, using summaryMetrics, + for self and child bundles. """ if self.summary_values is None: self.summary_values = {} if self.summary_metrics is not None: - # Build array of metric values, to use for (most) summary statistics. + # Build array of metric values, to use for summary statistics. for m in self.summary_metrics: summary_name = m.name summary_val = m.run(self.metric_values, self.slicer.slice_points["H"]) @@ -280,19 +285,20 @@ def __init__(self, bundle_dict, out_dir=".", results_db=None, verbose=True): def _check_compatible(self, metric_bundle1, metric_bundle2): """Check if two MetricBundles are "compatible". - Compatible indicates that the constraints, the slicers, and the maps are the same, and + Compatible indicates that the constraints, the slicers, + and the maps are the same, and that the stackers do not interfere with each other (i.e. are not trying to set the same column in different ways). Returns True if the MetricBundles are compatible, False if not. Parameters ---------- - metric_bundle1 : MetricBundle - metric_bundle2 : MetricBundle + metric_bundle1 : `MetricBundle` + metric_bundle2 : `MetricBundle` Returns ------- - bool + match : `bool` """ if metric_bundle1.constraint != metric_bundle2.constraint: return False @@ -302,24 +308,29 @@ def _check_compatible(self, metric_bundle1, metric_bundle2): return False for stacker in metric_bundle1.stacker_list: for stacker2 in metric_bundle2.stacker_list: - # If the stackers have different names, that's OK, and if they are identical, that's ok. + # If the stackers have different names, that's OK, + # and if they are identical, that's ok. if (stacker.__class__.__name__ == stacker2.__class__.__name__) & (stacker != stacker2): return False # But if we got this far, everything matches. return True def _find_compatible(self, test_keys): - """ "Private utility to find which metricBundles with keys in the list 'test_keys' can be calculated - at the same time -- having the same slicer, constraint, maps, and compatible stackers. + """ "Private utility to find which metricBundles with keys in the + list 'test_keys' can be calculated + at the same time -- having the same slicer, constraint, maps, + and compatible stackers. Parameters ----------- - test_keys : list - List of the dictionary keys (of self.bundle_dict) to test for compatibilility. + test_keys : `list` + List of the dictionary keys (of self.bundle_dict) to + test for compatibilility. Returns -------- list of lists - Returns test_keys, split into separate lists of compatible metricBundles. + Returns test_keys, split into separate lists of + compatible metricBundles. """ compatible_lists = [] for k in test_keys: @@ -334,14 +345,17 @@ def _find_compatible(self, test_keys): found_compatible = False checked_all = False while not (found_compatible) and not (checked_all): - # Go through the existing lists in compatible_lists, to see if this metricBundle matches. + # Go through the existing lists in compatible_lists, to see + # if this metricBundle matches. for compatible_list in compatible_lists: - # Compare to all the metricBundles in this subset, to check all stackers are compatible. + # Compare to all the metricBundles in this subset, + # to check all stackers are compatible. found_compatible = True for comparison_key in compatible_list: compatible = self._check_compatible(self.bundle_dict[comparison_key], b) if not compatible: - # Found a metricBundle which is not compatible, so stop and go onto the next subset. + # Found a metricBundle which is not compatible, + # so stop and go onto the next subset. found_compatible = False break checked_all = True @@ -356,14 +370,17 @@ def _find_compatible(self, test_keys): return compatible_lists def run_constraint(self, constraint): - """Calculate the metric values for all the metricBundles which match this constraint in the - metricBundleGroup. Also calculates child metrics and summary statistics, and writes all to disk. - (work is actually done in _runCompatible, so that only completely compatible sets of metricBundles + """Calculate the metric values for all the metricBundles which + match this constraint in the metricBundleGroup. + Also calculates child metrics and summary statistics, + and writes all to disk. + (work is actually done in _runCompatible, so that only completely + compatible sets of metricBundles run at the same time). Parameters ---------- - constraint : str + constraint : `str` SQL-where or pandas constraint for the metricBundles. """ # Find the dict keys of the bundles which match this constraint. @@ -376,7 +393,8 @@ def run_constraint(self, constraint): # Identify the observations which are relevant for this constraint. # This sets slicer.obs (valid for all H values). self.slicer.subset_obs(constraint) - # Identify the sets of these metricBundles can be run at the same time (also have the same stackers). + # Identify the sets of these metricBundles can be run at the same time + # (also have the same stackers). compatible_lists = self._find_compatible(keys_matching_constraint) # And now run each of those subsets of compatible metricBundles. @@ -384,22 +402,24 @@ def run_constraint(self, constraint): self._run_compatible(compatible_list) def _run_compatible(self, compatible_list): - """Calculate the metric values for set of (parent and child) bundles, as well as the summary stats, - and write to disk. + """Calculate the metric values for set of (parent and child) bundles, + as well as the summary stats, and write to disk. Parameters ----------- - compatible_list : list - List of dictionary keys, of the metricBundles which can be calculated together. - This means they are 'compatible' and have the same slicer, constraint, and non-conflicting - mappers and stackers. + compatible_list : `list` + List of dictionary keys, of the metricBundles which can be + calculated together. This means they are 'compatible' and have + the same slicer, constraint, and non-conflicting mappers and + stackers. """ if self.verbose: print("Running metrics %s" % compatible_list) - b_dict = self.bundle_dict # {key: self.bundle_dict.get(key) for key in compatible_list} + b_dict = self.bundle_dict - # Find the unique stackers and maps. These are already "compatible" (as id'd by compatible_list). + # Find the unique stackers and maps. + # These are already "compatible" (as id'd by compatible_list). uniq_stackers = [] all_stackers = [] uniq_maps = [] @@ -435,12 +455,14 @@ def _run_compatible(self, compatible_list): # Run all the parent metrics. for k in compatible_list: b = self.bundle_dict[k] - # Mask the parent metric (and then child metrics) if there was no data. + # Mask the parent metric (and then child metrics) + # if there was no data. if len(sso_obs) == 0: b.metric_values.mask[i][j] = True for cb in list(b.child_bundles.values()): cb.metric_values.mask[i][j] = True - # Otherwise, calculate the metric value for the parent, and then child. + # Otherwise, calculate the metric value for the parent, + # and then child. else: # Calculate for the parent. m_val = b.metric.run(sso_obs, slice_point["orbit"], Hval) @@ -449,7 +471,8 @@ def _run_compatible(self, compatible_list): b.metric_values.mask[i][j] = True for cb in b.child_bundles.values(): cb.metric_values.mask[i][j] = True - # Otherwise, set the parent value and calculate the child metric values as well. + # Otherwise, set the parent value and calculate + # the child metric values as well. else: b.metric_values.data[i][j] = m_val for cb in b.child_bundles.values(): @@ -469,9 +492,7 @@ def _run_compatible(self, compatible_list): b.write(out_dir=self.out_dir, results_db=self.results_db) def run_all(self): - """ - Run all constraints and metrics for these moMetricBundles. - """ + """Run all constraints and metrics for these moMetricBundles.""" for constraint in self.constraints: self.run_constraint(constraint) if self.verbose: @@ -487,7 +508,11 @@ def plot_all( closefigs=True, ): """ - Make a few generically desired plots. This needs more flexibility in the future. + Make a few generically desired plots. + Given the nature of the outputs for much of the moving object + metrics, a good deal of the plotting for the moving object batch + is handled in a custom manner joining together multiple + metricsbundles. """ plot_handler = PlotHandler( out_dir=self.out_dir, diff --git a/rubin_sim/maf/metrics/__init__.py b/rubin_sim/maf/metrics/__init__.py index 22c210eed..3c39eb153 100644 --- a/rubin_sim/maf/metrics/__init__.py +++ b/rubin_sim/maf/metrics/__init__.py @@ -17,7 +17,6 @@ from .hourglass_metric import * from .incremental_template_metric import * from .kuiper_metrics import * -from .long_gap_agn_metric import * from .mo_metrics import * from .mo_summary_metrics import * from .night_pointing_metric import * @@ -30,7 +29,6 @@ from .season_metrics import * from .simple_metrics import * from .sky_sat_metric import * -from .slew_metrics import * from .sn_cadence_metric import * from .sn_n_sn_metric import * from .sn_sl_metric import * diff --git a/rubin_sim/maf/metrics/agn_time_lag_metric.py b/rubin_sim/maf/metrics/agn_time_lag_metric.py index 12f368a37..b870f7dc3 100644 --- a/rubin_sim/maf/metrics/agn_time_lag_metric.py +++ b/rubin_sim/maf/metrics/agn_time_lag_metric.py @@ -50,7 +50,8 @@ def __init__( **kwargs, ) - # Calculate NQUIST value for time-lag and sampling time (redshift is included in formula if desired) + # Calculate NQUIST value for time-lag and sampling time + # (redshift is included in formula if desired) def _get_nquist_value(self, caden, lag, z): return lag / ((1 + z) * caden) diff --git a/rubin_sim/maf/metrics/agnstructure.py b/rubin_sim/maf/metrics/agnstructure.py index 6f54d5db2..3b95142c4 100644 --- a/rubin_sim/maf/metrics/agnstructure.py +++ b/rubin_sim/maf/metrics/agnstructure.py @@ -12,31 +12,35 @@ class SFUncertMetric(BaseMetric): - """Structure Function (SF) Uncertainty Metric. Developed on top of LogTGaps + """Structure Function (SF) Uncertainty Metric. + Developed on top of LogTGaps Adapted from Weixiang Yu & Gordon Richards at: - https://github.com/RichardsGroup/LSST_SF_Metric/blob/main/notebooks/00_SFErrorMetric.ipynb + https://github.com/RichardsGroup/ + LSST_SF_Metric/blob/main/notebooks/00_SFErrorMetric.ipynb Parameters ---------- - mag: `float` (22) + mag : `float` The magnitude of the fiducial object. Default 22. - times_col: `str` ('observationStartMJD') + times_col : `str` Time column name. Defaults to "observationStartMJD". - all_gaps: `bool` (True) + all_gaps : `bool` Whether to use all gaps (between any two pairs of observations). If False, only use consecutive paris. Defaults to True. - units: `str` ('mag') + units : `str` Unit of this metric. Defaults to "mag". - bins: `object` - An array of bin edges. Defaults to "np.logspace(0, np.log10(3650), 16)" for a + bins : `object` + An array of bin edges. + Defaults to "np.logspace(0, np.log10(3650), 16)" for a total of 15 (final) bins. - weight: `object + weight : `object` The weight assigned to each delta_t bin for deriving the final metric. - Defaults to flat weighting with sum of 1. Should have length 1 less than bins. - snr_cut : float (5) + Defaults to flat weighting with sum of 1. + Should have length 1 less than bins. + snr_cut : `float` Ignore observations below an SNR limit, default 5. - dust : `bool` (True) + dust : `bool` Apply dust extinction to the fiducial object magnitude. Default True. """ @@ -115,15 +119,17 @@ def run(self, data_slice, slice_point=None): else: dts = np.diff(times) - # bin delta_t using provided bins; if zero pair found at any delta_t bin, - # replace 0 with 0.01 to avoid the exploding 1/sqrt(n) term in this metric + # bin delta_t using provided bins; + # if zero pair found at any delta_t bin, + # replace 0 with 0.01 to avoid the exploding 1/sqrt(n) term + # in this metric result, bins = np.histogram(dts, self.bins) new_result = np.where(result > 0, result, 0.01) # compute photometric_error^2 population variance and population mean # note that variance is replaced by median_absolute_deviate^2 - # mean is replaced by median in this implementation to make it robust to - # outliers in simulations (e.g., dcr simulations) + # mean is replaced by median in this implementation to make it robust + # to outliers in simulations (e.g., dcr simulations) err_var = mag_err**2 err_var_mu = np.median(err_var) err_var_std = mad_std(err_var) diff --git a/rubin_sim/maf/metrics/area_summary_metrics.py b/rubin_sim/maf/metrics/area_summary_metrics.py index 78119965d..1890953c8 100644 --- a/rubin_sim/maf/metrics/area_summary_metrics.py +++ b/rubin_sim/maf/metrics/area_summary_metrics.py @@ -8,21 +8,23 @@ class AreaSummaryMetric(BaseMetric): """ - Find the min/max of a value in the best area. This is a handy substitute for when - users want to know "the WFD value". + Find the min/max of a value over the area with the 'best' results + in the metric. + This is a handy substitute for when users want to know "the WFD value". Parameters ---------- - area : float (18000) + area : `float` The area to consider (sq degrees) - decreasing : bool (True) - Should the values be sorted by increasing or decreasing order. For values where - "larger is better", decreasing is probably what you want. For metrics where - "smaller is better" (e.g., astrometric precission), set decreasing to False. + decreasing : `bool` + Should the values be sorted by increasing or decreasing order. + For values where "larger is better", decreasing (True) is probably + what you want. For metrics where "smaller is better" + (e.g., astrometric precission), set decreasing to False. reduce_func : None - The function to reduce the clipped values by. Will default to min/max depending on - the bool val of the decreasing kwarg. - + The function to reduce the clipped values by. + Will default to min/max depending on the bool val of the decreasing + kwarg. """ def __init__( @@ -64,11 +66,12 @@ def run(self, data_slice, slice_point=None): class AreaThresholdMetric(BaseMetric): - """ - Find the amount of area on the sky that meets a given threshold value. + """Find the amount of area on the sky that meets a given threshold value. - The area per pixel is determined from the size of the metric_values array passed to the summary metric. - This assumes that both all values are passed and that the metric was calculated with a healpix slicer. + The area per pixel is determined from the size of the metric_values + array passed to the summary metric. + This assumes that both all values are passed and that the metric was + calculated with a healpix slicer. Parameters ---------- diff --git a/rubin_sim/maf/metrics/base_metric.py b/rubin_sim/maf/metrics/base_metric.py index bb6fb2b63..030c07aba 100644 --- a/rubin_sim/maf/metrics/base_metric.py +++ b/rubin_sim/maf/metrics/base_metric.py @@ -1,10 +1,4 @@ # Base class for metrics - defines methods which must be implemented. -# If a metric calculates a vector or list at each gridpoint, then there -# should be additional 'reduce_*' functions defined, to convert the vector -# into scalar (and thus plottable) values at each gridpoint. -# The philosophy behind keeping the vector instead of the scalar at each gridpoint -# is that these vectors may be expensive to compute; by keeping/writing the full -# vector we permit multiple 'reduce' functions to be executed on the same data. __all__ = ("MetricRegistry", "BaseMetric", "ColRegistry") @@ -114,8 +108,9 @@ class BaseMetric(metaclass=MetricRegistry): ---------- col : `str` or `list` [`str`] Names of the data columns that the metric will use. - The columns required for each metric is tracked in the ColRegistry, and used to retrieve data - from the opsim database. Can be a single string or a list. + The columns required for each metric is tracked in the ColRegistry, + and used to retrieve data from the opsim database. + Can be a single string or a list. metric_name : `str` Name to use for the metric (optional - if not set, will be derived). maps : `list` [`rubin_sim.maf.maps`] @@ -143,9 +138,11 @@ def __init__( badval=-666, mask_val=None, ): - # Turn cols into numpy array so we know we can iterate over the columns. + # Turn cols into numpy array so we know + # we can iterate over the columns. self.col_name_arr = np.array(col, copy=False, ndmin=1) - # To support simple metrics operating on a single column, set self.colname + # To support simple metrics operating on a single column, + # set self.colname if len(self.col_name_arr) == 1: self.colname = self.col_name_arr[0] # Add the columns to the colRegistry. @@ -161,7 +158,8 @@ def __init__( # Save a unique name for the metric. self.name = metric_name if self.name is None: - # If none provided, construct our own from the class name and the data columns. + # If none provided, construct our own from the class name + # and the data columns. self.name = ( self.__class__.__name__.replace("Metric", "", 1) + " " @@ -197,7 +195,7 @@ def __init__( # Default to only return one metric value per slice self.shape = 1 - def run(self, data_slice, slice_point=None): + def run(self, data_slice, slice_point): """Calculate metric values. Parameters @@ -207,11 +205,11 @@ def run(self, data_slice, slice_point=None): use to calculate metric values at each slice_point. slice_point : `dict` or None Dictionary of slice_point metadata passed to each metric. - E.g. the ra/dec of the healpix pixel or opsim fieldId. + E.g. the ra/dec of the healpix pixel. Returns ------- - metricValue: `int` `float` or `object` + metricValue : `int` `float` or `object` The metric value at each slice_point. """ raise NotImplementedError("Please implement your metric calculation.") diff --git a/rubin_sim/maf/metrics/brown_dwarf_metric.py b/rubin_sim/maf/metrics/brown_dwarf_metric.py index b652466c0..6fc299858 100644 --- a/rubin_sim/maf/metrics/brown_dwarf_metric.py +++ b/rubin_sim/maf/metrics/brown_dwarf_metric.py @@ -43,32 +43,42 @@ def bd_colors(spec_type): class BDParallaxMetric(BaseMetric): - """Calculate the distance to which one could reach a parallax SNR for a given object + """Calculate the distance to which one could reach a parallax SNR for a + given object + Modification of ParallaxMetric, illustrated in - https://github.com/jgizis/LSST-BD-Cadence/blob/main/bd_allLT_baseline_17.ipynb + https://github.com/jgizis/ + LSST-BD-Cadence/blob/main/bd_allLT_baseline_17.ipynb - Uses columns ra_pi_amp and dec_pi_amp, calculated by the ParallaxFactorStacker. + Uses columns ra_pi_amp and dec_pi_amp, + calculated by the ParallaxFactorStacker. Parameters ---------- metricName : `str`, opt Default 'parallax'. m5_col : `str`, opt - The default column name for m5 information in the input data. Default fiveSigmaDepth. + The default column name for m5 information in the input data. + Default fiveSigmaDepth. filter_col : `str`, opt The column name for the filter information. Default filter. seeing_col : `str`, opt - The column name for the seeing information. Since the astrometry errors are based on the physical - size of the PSF, this should be the FWHM of the physical psf. Default seeingFwhmGeom. - mags : `dict` (None) - The absolute magnitude of the obeject in question. Keys of filter name, values in mags. + The column name for the seeing information. + Since the astrometry errors are based on the physical + size of the PSF, this should be the FWHM of the physical psf. + Default seeingFwhmGeom. + mags : `dict` or None + The absolute magnitude of the object in question. + Keys of filter name, values in mags. Defaults to an L7 spectral type if None. - distances : np.array + distances : `np.array`, (N,) Distances to try putting the object at (pc). atm_err : `float`, opt - The expected centroiding error due to the atmosphere, in arcseconds. Default 0.01. + The expected centroiding error due to the atmosphere, in arcseconds. + Default 0.01. badval : `float`, opt - The value to return when the metric value cannot be calculated. Default 0. + The value to return when the metric value cannot be calculated. + Default 0. """ def __init__( @@ -107,7 +117,8 @@ def __init__( def _final_sigma(self, position_errors, ra_pi_amp, dec_pi_amp): """Assume parallax in RA and DEC are fit independently, then combined. - All inputs assumed to be arcsec""" + All inputs assumed to be arcsec. + """ sigma_a = position_errors / ra_pi_amp sigma_b = position_errors / dec_pi_amp sigma_ra = np.sqrt(1.0 / np.sum(1.0 / sigma_a**2, axis=1)) @@ -139,7 +150,7 @@ def run(self, data_slice, slice_point=None): class VolumeSumMetric(BaseMetric): - """Compute the total volume assuming a metric has values of distance""" + """Compute the total volume assuming a metric has values of distance.""" def __init__(self, col=None, metric_name="VolumeSum", nside=None, **kwargs): super(VolumeSumMetric, self).__init__(col=col, metric_name=metric_name, **kwargs) diff --git a/rubin_sim/maf/metrics/cadence_metrics.py b/rubin_sim/maf/metrics/cadence_metrics.py index 85c5b8356..385d22a46 100644 --- a/rubin_sim/maf/metrics/cadence_metrics.py +++ b/rubin_sim/maf/metrics/cadence_metrics.py @@ -24,14 +24,17 @@ def __init__(self, filter_col="filter", metric_name="fS", **kwargs): super().__init__(cols=cols, metric_name=metric_name, units="fS", **kwargs) def run(self, data_slice, slice_point=None): - # We could import this from the m5_flat_sed values, but it makes sense to calculate the m5 - # directly from the throughputs. This is easy enough to do and will allow variation of + # We could import this from the m5_flat_sed values, + # but it makes sense to calculate the m5 + # directly from the throughputs. This is easy enough to do and + # will allow variation of # the throughput curves and readnoise and visit length, etc. pass class TemplateExistsMetric(BaseMetric): - """Calculate the fraction of images with a previous template image of desired quality.""" + """Calculate the fraction of images with a previous template + image of desired quality.""" def __init__( self, @@ -52,7 +55,8 @@ def run(self, data_slice, slice_point=None): data_slice.sort(order=self.observation_start_mjd_col) # Find the minimum seeing up to a given time seeing_mins = np.minimum.accumulate(data_slice[self.seeing_col]) - # Find the difference between the seeing and the minimum seeing at the previous visit + # Find the difference between the seeing and the minimum seeing + # at the previous visit seeing_diff = data_slice[self.seeing_col] - np.roll(seeing_mins, 1) # First image never has a template; check how many others do good = np.where(seeing_diff[1:] >= 0.0)[0] @@ -63,14 +67,16 @@ def run(self, data_slice, slice_point=None): class UniformityMetric(BaseMetric): """Calculate how uniformly the observations are spaced in time. - This is based on how a KS-test works: look at the cumulative distribution of observation dates, + This is based on how a KS-test works: + look at the cumulative distribution of observation dates, and compare to a perfectly uniform cumulative distribution. Perfectly uniform observations = 0, perfectly non-uniform = 1. Parameters ---------- mjd_col : `str`, optional - The column containing time for each observation. Default "observationStartMJD". + The column containing time for each observation. + Default "observationStartMJD". survey_length : `float`, optional The overall duration of the survey. Default 10. """ @@ -85,7 +91,8 @@ def run(self, data_slice, slice_point=None): # If only one observation, there is no uniformity if data_slice[self.mjd_col].size == 1: return 1 - # Scale dates to lie between 0 and 1, where 0 is the first observation date and 1 is surveyLength + # Scale dates to lie between 0 and 1, + # where 0 is the first observation date and 1 is surveyLength dates = (data_slice[self.mjd_col] - data_slice[self.mjd_col].min()) / (self.survey_length * 365.25) dates.sort() # Just to be sure n_cum = np.arange(1, dates.size + 1) / float(dates.size) @@ -96,7 +103,8 @@ def run(self, data_slice, slice_point=None): class GeneralUniformityMetric(BaseMetric): """Calculate how uniformly any values are distributed. - This is based on how a KS-test works: look at the cumulative distribution of data, + This is based on how a KS-test works: + look at the cumulative distribution of data, and compare to a perfectly uniform cumulative distribution. Perfectly uniform observations = 0, perfectly non-uniform = 1. To be "perfectly uniform" here, the endpoints need to be included. @@ -105,7 +113,8 @@ class GeneralUniformityMetric(BaseMetric): ---------- col : `str`, optional The column of data to use for the metric. - The default is "observationStartMJD" as this is most typically used with time. + The default is "observationStartMJD" as this is most + typically used with time. min_value : `float`, optional The minimum value expected for the data. Default None will calculate use the minimum value in this dataslice @@ -117,7 +126,6 @@ class GeneralUniformityMetric(BaseMetric): """ def __init__(self, col="observationStartMJD", units="", min_value=None, max_value=None, **kwargs): - """survey_length = time span of survey (years)""" self.col = col super().__init__(col=self.col, units=units, **kwargs) self.min_value = min_value @@ -127,7 +135,8 @@ def run(self, data_slice, slice_point=None): # If only one observation, there is no uniformity if data_slice[self.col].size == 1: return 1 - # Scale values to lie between 0 and 1, where 0 is the min_value and 1 is max_value + # Scale values to lie between 0 and 1, + # where 0 is the min_value and 1 is max_value if self.min_value is None: min_value = data_slice[self.col].min() else: @@ -144,9 +153,11 @@ def run(self, data_slice, slice_point=None): class RapidRevisitUniformityMetric(BaseMetric): - """Calculate uniformity of time between consecutive visits on short timescales (for RAV1). + """Calculate uniformity of time between consecutive visits on + short timescales (for RAV1). - Uses a the same 'uniformity' calculation as the UniformityMetric, based on the KS-test. + Uses the same 'uniformity' calculation as the UniformityMetric, + based on the KS-test. A value of 0 is perfectly uniform; a value of 1 is purely non-uniform. Parameters @@ -154,7 +165,8 @@ class RapidRevisitUniformityMetric(BaseMetric): mjd_col : `str`, optional The column containing the 'time' value. Default observationStartMJD. min_nvisits : `int`, optional - The minimum number of visits required within the time interval (d_tmin to d_tmax). + The minimum number of visits required within the + time interval (d_tmin to d_tmax). Default 100. d_tmin : `float`, optional The minimum dTime to consider (in days). Default 40 seconds. @@ -176,7 +188,8 @@ def __init__( self.d_tmin = d_tmin self.d_tmax = d_tmax super().__init__(col=self.mjd_col, metric_name=metric_name, **kwargs) - # Update min_nvisits, as 0 visits will crash algorithm and 1 is nonuniform by definition. + # Update min_nvisits, as 0 visits will crash algorithm + # and 1 is nonuniform by definition. if self.min_nvisits <= 1: self.min_nvisits = 2 @@ -230,16 +243,19 @@ def run(self, data_slice, slice_point=None): class NRevisitsMetric(BaseMetric): - """Calculate the number of consecutive visits with time differences less than d_t. + """Calculate the number of consecutive visits with + time differences less than d_t. Parameters ---------- d_t : `float`, optional The time interval to consider (in minutes). Default 30. normed : `bool`, optional - Flag to indicate whether to return the total number of consecutive visits with time - differences less than d_t (False), or the fraction of overall visits (True). - Note that we would expect (if all visits occur in pairs within d_t) this fraction would be 0.5! + Flag to indicate whether to return the total number of + consecutive visits with time differences less than d_t (False), + or the fraction of overall visits (True). + Note that we would expect (if all visits occur in pairs within d_t) + this fraction would be 0.5! """ def __init__(self, mjd_col="observationStartMJD", d_t=30.0, normed=False, metric_name=None, **kwargs): @@ -267,12 +283,14 @@ def run(self, data_slice, slice_point=None): class IntraNightGapsMetric(BaseMetric): """ - Calculate the (reduce_func) of the gap between consecutive observations within a night, in hours. + Calculate the (reduce_func) of the gap between consecutive + observations within a night, in hours. Parameters ---------- reduce_func : function, optional - Function that can operate on array-like structures. Typically numpy function. + Function that can operate on array-like structures. + Typically numpy function. Default np.median. """ @@ -306,12 +324,14 @@ def run(self, data_slice, slice_point=None): class InterNightGapsMetric(BaseMetric): - """Calculate the (reduce_func) of the gap between consecutive observations in different nights, in days. + """Calculate the (reduce_func) of the gap between consecutive + observations in different nights, in days. Parameters ---------- reduce_func : function, optional - Function that can operate on array-like structures. Typically numpy function. + Function that can operate on array-like structures. + Typically numpy function. Default np.median. """ @@ -344,16 +364,18 @@ def run(self, data_slice, slice_point=None): class VisitGapMetric(BaseMetric): - """Calculate the (reduce_func) of the gap between any consecutive observations, in hours, - regardless of night boundaries. + """Calculate the (reduce_func) of the gap between any + consecutive observations, in hours, regardless of night boundaries. - Different from inter-night and intra-night gaps, between this is really just counting - all of the times between consecutive observations (not time between nights or time within a night). + Different from inter-night and intra-night gaps, + because this is really just counting all of the times between consecutive + observations (not time between nights or time within a night). Parameters ---------- reduce_func : function, optional - Function that can operate on array-like structures. Typically numpy function. + Function that can operate on array-like structures. + Typically numpy function. Default np.median. """ diff --git a/rubin_sim/maf/metrics/calibration_metrics.py b/rubin_sim/maf/metrics/calibration_metrics.py index 2c7484808..498f447de 100644 --- a/rubin_sim/maf/metrics/calibration_metrics.py +++ b/rubin_sim/maf/metrics/calibration_metrics.py @@ -97,11 +97,11 @@ def __init__( "This normalized version of the metric displays the " "estimated uncertainty in the parallax measurement, " ) - self.comment += "divided by the minimum parallax uncertainty possible " "(if all visits were six " - self.comment += ( - "months apart). Values closer to 1 indicate more optimal " - "scheduling for parallax measurement." - ) + self.comment += "divided by the minimum parallax uncertainty possible " + self.comment += "(if all visits were six months apart). " + self.comment += "Values closer to 1 indicate more optimal "\ + "scheduling for parallax measurement." + def _final_sigma(self, position_errors, ra_pi_amp, dec_pi_amp): """Assume parallax in RA and DEC are fit independently, then combined. @@ -127,7 +127,8 @@ def run(self, data_slice, slice_point=None): position_errors = mafUtils.astrom_precision(data_slice[self.seeing_col], snr, self.atm_err) sigma = self._final_sigma(position_errors, data_slice["ra_pi_amp"], data_slice["dec_pi_amp"]) if self.normalize: - # Leave the dec parallax as zero since one can't have ra and dec maximized at the same time. + # Leave the dec parallax as zero since one can't have + # ra and dec maximized at the same time. sigma = ( self._final_sigma( position_errors, @@ -146,32 +147,41 @@ class ProperMotionMetric(BaseMetric): Parameters ---------- - metricName : str, optional + metricName : `str`, optional Default 'properMotion'. - m5_col : str, optional - The default column name for m5 information in the input data. Default fiveSigmaDepth. - mjd_col : str, optional + m5_col : `str`, optional + The default column name for m5 information in the input data. + Default fiveSigmaDepth. + mjd_col : `str`, optional The column name for the exposure time. Default observationStartMJD. - filterCol : str, optional + filterCol : `str`, optional The column name for the filter information. Default filter. - seeing_col : str, optional - The column name for the seeing information. Since the astrometry errors are based on the physical - size of the PSF, this should be the FWHM of the physical psf. Default seeingFwhmGeom. - rmag : float, optional - The r magnitude of the fiducial star in r band. Other filters are sclaed using sedTemplate keyword. + seeing_col : `str`, optional + The column name for the seeing information. + Since the astrometry errors are based on the physical + size of the PSF, this should be the FWHM of the physical psf. + Default seeingFwhmGeom. + rmag : `float`, optional + The r magnitude of the fiducial star in r band. + Other filters are sclaed using sedTemplate keyword. Default 20.0 - SedTemplate : str, optional - The template to use. This can be 'flat' or 'O','B','A','F','G','K','M'. Default flat. - atm_err : float, optional - The expected centroiding error due to the atmosphere, in arcseconds. Default 0.01. + SedTemplate : `str`, optional + The template to use. This can be 'flat' or 'O','B','A','F','G','K','M'. + Default flat. + atm_err : `float`, optional + The expected centroiding error due to the atmosphere, in arcseconds. + Default 0.01. normalize : `bool`, optional - Compare the astrometric uncertainty to the uncertainty that would result if half the observations - were taken at the start and half at the end. A perfect survey will have a value close to 1, while + Compare the astrometric uncertainty to the uncertainty that would + result if half the observations were taken at the start and half + at the end. A perfect survey will have a value close to 1, while a poorly scheduled survey will be close to 0. Default False. - baseline : float, optional - The length of the survey used for the normalization, in years. Default 10. - badval : float, optional - The value to return when the metric value cannot be calculated. Default -666. + baseline : `float`, optional + The length of the survey used for the normalization, in years. + Default 10. + badval : `float`, optional + The value to return when the metric value cannot be calculated. + Default -666. """ def __init__( @@ -225,8 +235,8 @@ def __init__( self.comment += ( "This normalized version of the metric represents " "the estimated uncertainty in the proper " ) - self.comment += "motion divided by the minimum uncertainty possible " "(if all visits were " - self.comment += "obtained on the first and last days of the survey). " + self.comment += "motion divided by the minimum uncertainty possible " + self.comment += "(if all visits were obtained on the first and last days of the survey). " self.comment += "Values closer to 1 indicate more optimal scheduling." def run(self, data_slice, slice_point=None): @@ -257,50 +267,61 @@ def run(self, data_slice, slice_point=None): class ParallaxCoverageMetric(BaseMetric): - """ - Check how well the parallax factor is distributed. Subtracts the weighted mean position of the + """Check how well the parallax factor is distributed. + + Subtracts the weighted mean position of the parallax offsets, then computes the weighted mean radius of the points. - If points are well distributed, the mean radius will be near 1. If phase coverage is bad, - radius will be close to zero. + If points are well distributed, the mean radius will be near 1. + If phase coverage is bad, radius will be close to zero. - For points on the Ecliptic, uniform sampling should result in a metric value of ~0.5. + For points on the Ecliptic, uniform sampling should result in a + metric value of ~0.5. At the poles, uniform sampling would result in a metric value of ~1. - Conceptually, it is helpful to remember that the parallax motion of a star at the pole is - a (nearly circular) ellipse while the motion of a star on the ecliptic is a straight line. Thus, any - pair of observations separated by 6 months will give the full parallax range for a star on the pole - but only observations on very specific dates will give the full range for a star on the ecliptic. + Conceptually, it is helpful to remember that the parallax motion of a + star at the pole is a (nearly circular) ellipse while the motion of a + star on the ecliptic is a straight line. Thus, any pair of observations + separated by 6 months will give the full parallax range for a star on + the pole but only observations on very specific dates will give the + full range for a star on the ecliptic. - Optionally also demand that there are observations above the snr_limit kwarg spanning theta_range radians. + Optionally also demand that there are observations above the snr_limit + kwarg spanning theta_range radians. Parameters ---------- - m5_col: str, optional + m5_col : `str`, optional Column name for individual visit m5. Default fiveSigmaDepth. - mjd_col: str, optional + mjd_col : `str`, optional Column name for exposure time dates. Default observationStartMJD. - filter_col: str, optional + filter_col : `str`, optional Column name for filter. Default filter. - seeing_col: str, optional + seeing_col : `str`, optional Column name for seeing (assumed FWHM). Default seeingFwhmGeom. - rmag: float, optional - Magnitude of fiducial star in r filter. Other filters are scaled using sedTemplate keyword. + rmag : `float`, optional + Magnitude of fiducial star in r filter. + Other filters are scaled using sedTemplate keyword. Default 20.0 - sedTemplate: str, optional - Template to use (can be 'flat' or 'O','B','A','F','G','K','M'). Default 'flat'. - atm_err: float, optional - Centroiding error due to atmosphere in arcsec. Default 0.01 (arcseconds). - theta_range: float, optional - Range of parallax offset angles to demand (in radians). Default=0 (means no range requirement). - snr_limit: float, optional - Only include points above the snr_limit when computing theta_range. Default 5. + sedTemplate : `str`, optional + Template to use (can be 'flat' or 'O','B','A','F','G','K','M'). + Default 'flat'. + atm_err : `float`, optional + Centroiding error due to atmosphere in arcsec. + Default 0.01 (arcseconds). + theta_range : `float`, optional + Range of parallax offset angles to demand (in radians). + Default=0 (means no range requirement). + snr_limit : `float`, optional + Only include points above the snr_limit when computing theta_range. + Default 5. Returns -------- - metricValue: float + metricValu e: `float` Returns a weighted mean of the length of the parallax factor vectors. Values near 1 imply that the points are well distributed. Values near 0 imply that the parallax phase coverage is bad. - Near the ecliptic, uniform sampling results in metric values of about 0.5. + Near the ecliptic, uniform sampling results in metric values + of about 0.5. Notes ----- @@ -354,7 +375,7 @@ def _theta_check(self, ra_pi_amp, dec_pi_amp, snr): theta = theta - np.min(theta) result = 0.0 if np.max(theta) >= self.theta_range: - # Check that things are in differnet quadrants + # Check that things are in different quadrants theta = (theta + np.pi) % 2.0 * np.pi theta = theta - np.min(theta) if np.max(theta) >= self.theta_range: @@ -397,36 +418,42 @@ def run(self, data_slice, slice_point=None): class ParallaxDcrDegenMetric(BaseMetric): - """Use the full parallax and DCR displacement vectors to find if they are degenerate. + """Use the full parallax and DCR displacement vectors to find if they + are degenerate. Parameters ---------- - metricName : str, optional + metricName : `str`, optional Default 'ParallaxDcrDegenMetric'. - seeing_col : str, optional + seeing_col : `str`, optional Default 'FWHMgeom' - m5_col : str, optional + m5_col : `str`, optional Default 'fiveSigmaDepth' - filter_col : str + filter_col : `str` Default 'filter' - atm_err : float - Minimum error in photometry centroids introduced by the atmosphere (arcseconds). Default 0.01. - rmag : float + atm_err : `float` + Minimum error in photometry centroids introduced by the atmosphere + (arcseconds). Default 0.01. + rmag : `float` r-band magnitude of the fiducual star that is being used (mag). - SedTemplate : str - The SED template to use for fiducia star colors, passed to rubin_scheduler.utils.stellarMags. + SedTemplate : `str` + The SED template to use for fiducia star colors, + passed to rubin_scheduler.utils.stellarMags. Default 'flat' - tol : float - Tolerance for how well curve_fit needs to work before believing the covariance result. + tol : `float` + Tolerance for how well curve_fit needs to work before + believing the covariance result. Default 0.05. Returns ------- - metricValue : float - Returns the correlation coefficient between the best-fit parallax amplitude and DCR amplitude. - The RA and Dec offsets are fit simultaneously. Values close to zero are good, values close to +/- 1 - are bad. Experience with fitting Monte Carlo simulations suggests the astrometric fits start - becoming poor around a correlation of 0.7. + metricValue : `float` + Returns the correlation coefficient between the best-fit parallax + amplitude and DCR amplitude. + The RA and Dec offsets are fit simultaneously. + Values close to zero are good, values close to +/- 1 are bad. + Experience with fitting Monte Carlo simulations suggests the + astrometric fits start becoming poor around a correlation of 0.7. """ def __init__( @@ -465,39 +492,46 @@ def __init__( self.atm_err = atm_err def _positions(self, x, a, b): - """ - Function to find parallax and dcr amplitudes + """Function to find parallax and dcr amplitudes - x should be a vector with [[parallax_x1, parallax_x2..., parallax_y1, parallax_y2...], + x should be a vector with [[parallax_x1, parallax_x2..., + parallax_y1, parallax_y2...], [dcr_x1, dcr_x2..., dcr_y1, dcr_y2...]] """ result = a * x[0, :] + b * x[1, :] return result def run(self, data_slice, slice_point=None): - # The idea here is that we calculate position errors (in RA and Dec) for all observations. - # Then we generate arrays of the parallax offsets (delta RA parallax = ra_pi_amp, etc) - # and the DCR offsets (delta RA DCR = ra_dcr_amp, etc), and just add them together into one - # RA (and Dec) offset. Then, we try to fit for how we combined these offsets, but while - # considering the astrometric noise. If we can figure out that we just added them together - # (i.e. the curve_fit result is [a=1, b=1] for the function _positions above) - # then we should be able to disentangle the parallax and DCR offsets when fitting 'for real'. + # The idea here is that we calculate position errors (in RA and Dec) + # for all observations. Then we generate arrays of the parallax + # offsets (delta RA parallax = ra_pi_amp, etc) and the DCR offsets + # (delta RA DCR = ra_dcr_amp, etc), and just add them together into one + # RA (and Dec) offset. Then, we try to fit for how we combined these + # offsets, but while considering the astrometric noise. If we can + # figure out that we just added them together + # (i.e. the curve_fit result is [a=1, b=1] for the function + # _positions above) then we should be able to disentangle the + # parallax and DCR offsets when fitting 'for real'. # compute SNR for all observations snr = np.zeros(len(data_slice), dtype="float") for filt in np.unique(data_slice[self.filter_col]): in_filt = np.where(data_slice[self.filter_col] == filt) snr[in_filt] = mafUtils.m52snr(self.mags[filt], data_slice[self.m5_col][in_filt]) # Compute the centroiding uncertainties - # Note that these centroiding uncertainties depend on the physical size of the PSF, thus - # we are using seeingFwhmGeom for these metrics, not seeingFwhmEff. + # Note that these centroiding uncertainties depend on the physical + # size of the PSF, thus we are using seeingFwhmGeom for these metrics, + # not seeingFwhmEff. position_errors = mafUtils.astrom_precision(data_slice[self.seeing_col], snr, self.atm_err) - # Construct the vectors of RA/Dec offsets. xdata is the "input data". ydata is the "output". + # Construct the vectors of RA/Dec offsets. xdata is the "input data". + # ydata is the "output". xdata = np.empty((2, data_slice.size * 2), dtype=float) xdata[0, :] = np.concatenate((data_slice["ra_pi_amp"], data_slice["dec_pi_amp"])) xdata[1, :] = np.concatenate((data_slice["ra_dcr_amp"], data_slice["dec_dcr_amp"])) ydata = np.sum(xdata, axis=0) - # Use curve_fit to compute covariance between parallax and dcr amplitudes - # Set the initial guess slightly off from the correct [1,1] to make sure it iterates. + # Use curve_fit to compute covariance between parallax + # and dcr amplitudes + # Set the initial guess slightly off from the correct [1,1] to + # make sure it iterates. popt, pcov = curve_fit( self._positions, xdata, @@ -511,7 +545,8 @@ def run(self, data_slice, slice_point=None): return self.badval # Covariance between best fit parallax amplitude and DCR amplitude. cov = pcov[1, 0] - # Convert covarience between parallax and DCR amplitudes to normalized correlation + # Convert covariance between parallax and DCR amplitudes to normalized + # correlation perr = np.sqrt(np.diag(pcov)) correlation = cov / (perr[0] * perr[1]) result = correlation @@ -522,13 +557,26 @@ def run(self, data_slice, slice_point=None): def calc_dist_cosines(ra1, dec1, ra2, dec2): - # Taken from simSelfCalib.py """Calculates distance on a sphere using spherical law of cosines. - Give this function RA/Dec values in radians. Returns angular distance(s), in radians. - Note that since this is all numpy, you could input arrays of RA/Decs.""" + Note: floats can be replaced by numpy arrays of RA/Dec. + For very small distances, rounding errors may cause distance errors. + + Parameters + ---------- + ra1, dec1 : `float`, `float` + RA and Dec of one point. (radians) + ra2, dec2 : `float`, `float` + RA and Dec of another point. (radians) + + Returns + ------- + distance : `float` + Angular distance between the points in radians. + """ # This formula can have rounding errors for case where distances are small. - # Oh, the joys of wikipedia - http://en.wikipedia.org/wiki/Great-circle_distance + # Oh, the joys of wikipedia - + # http://en.wikipedia.org/wiki/Great-circle_distance # For the purposes of these calculations, this is probably accurate enough. D = np.sin(dec2) * np.sin(dec1) + np.cos(dec1) * np.cos(dec2) * np.cos(ra2 - ra1) D = np.arccos(D) @@ -536,7 +584,9 @@ def calc_dist_cosines(ra1, dec1, ra2, dec2): class RadiusObsMetric(BaseMetric): - """find the radius in the focal plane. returns things in degrees.""" + """Evaluate slice point radial position in the focal plane of each visit, + reducing to the mean, rms and full range of these radial distances. + """ def __init__( self, metric_name="radiusObs", ra_col="fieldRA", dec_col="fieldDec", units="radians", **kwargs diff --git a/rubin_sim/maf/metrics/coverage_metric.py b/rubin_sim/maf/metrics/coverage_metric.py index 23ef98e7e..d20333680 100644 --- a/rubin_sim/maf/metrics/coverage_metric.py +++ b/rubin_sim/maf/metrics/coverage_metric.py @@ -6,22 +6,26 @@ class YearCoverageMetric(BaseMetric): - """Count the number of bins covered by night_col -- default bins are 'years'. - Handy for checking that a point on the sky gets observed every year, as the default settings - result in the metric returning the number years in the data_slice (when used with a HealpixSlicer). + """Count the number of `bins` covered by night_col. + + The default `bins` cover years 0 to 10. + Handy for checking that a point on the sky gets observed every year, + as the default settings result in the metric returning the number years + in the data_slice (when used with a HealpixSlicer). Parameters ---------- - night_col: str, optional + night_col : `str`, opt Data column to histogram. Default 'night'. - bins: numpy.ndarray, optional - Bins to use in the histogram. Default corresponds to years 0-10 (with 365.25 nights per year). - units: str, optional + bins : `np.ndarray`, (N,), opt + Bins to use in the histogram. Default corresponds to years 0-10 + (with 365.25 nights per year). + units : `str`, opt Units to use for the metric result. Default 'N years'. Returns ------- - integer + nbins : `int` Number of histogram bins where the histogram value is greater than 0. Typically this will be the number of years in the 'night_col'. """ diff --git a/rubin_sim/maf/metrics/crowding_metric.py b/rubin_sim/maf/metrics/crowding_metric.py index ddb9e1fdf..4e3e7caff 100644 --- a/rubin_sim/maf/metrics/crowding_metric.py +++ b/rubin_sim/maf/metrics/crowding_metric.py @@ -6,33 +6,35 @@ from rubin_sim.maf.metrics import BaseMetric -# Modifying from Knut Olson's fork at: -# https://github.com/knutago/sims_maf_contrib/blob/master/tutorials/CrowdingMetric.ipynb +# Originally contributed by Knut Olson (@knutago). def _comp_crowd_error(mag_vector, lum_func, seeing, single_mag=None): """ - Compute the photometric crowding error given the luminosity function and best seeing. + Compute the photometric crowding error given the luminosity + function and best seeing. + + Equation from Olsen, Blum, & Rigaut 2003, AJ, 126, 452 Parameters ---------- - mag_vector : np.array + mag_vector : `np.array` (N,) Stellar magnitudes. - lum_func : np.array + lum_func : `np.array` (N,) Stellar luminosity function. - seeing : float - The best seeing conditions. Assuming forced-photometry can use the best seeing conditions + seeing : `float` + The best seeing conditions. + Assuming forced-photometry can use the best seeing conditions to help with confusion errors. - single_mag : float (None) - If single_mag is None, the crowding error is calculated for each mag in mag_vector. If - single_mag is a float, the crowding error is interpolated to that single value. + single_mag : `float` or None + If single_mag is None, the crowding error is calculated + for each mag in mag_vector. If single_mag is a float, + the crowding error is interpolated to that single value. Returns ------- - np.array + mag_uncert : `np.array` (N,) Magnitude uncertainties. - - Equation from Olsen, Blum, & Rigaut 2003, AJ, 126, 452 """ lum_area_arcsec = 3600.0**2 lum_vector = 10 ** (-0.4 * mag_vector) @@ -47,7 +49,29 @@ def _comp_crowd_error(mag_vector, lum_func, seeing, single_mag=None): class CrowdingM5Metric(BaseMetric): - """Return the magnitude at which the photometric error exceeds crowding_error threshold.""" + """Calculate the magnitude at which the photometric error exceeds + the crowding error threshold. + + Parameters + ---------- + crowding_error : `float`, optional + The magnitude uncertainty from crowding in magnitudes. + Default 0.1 mags. + filtername : `str`, optional + The bandpass in which to calculate the crowding limit. Default r. + seeing_col : `str`, optional + The name of the seeing column. + m5Col : `str`, optional + The name of the m5 depth column. + maps : `list` [`str`], optional + Names of maps required for the metric. + + Returns + ------- + mag : `float` + The magnitude of a star which has a photometric error of + `crowding_error` + """ def __init__( self, @@ -58,26 +82,6 @@ def __init__( maps=["StellarDensityMap"], **kwargs, ): - """ - Parameters - ---------- - crowding_error : float, optional - The magnitude uncertainty from crowding in magnitudes. Default 0.1 mags. - filtername: str, optional - The bandpass in which to calculate the crowding limit. Default r. - seeing_col : str, optional - The name of the seeing column. - m5Col : str, optional - The name of the m5 depth column. - maps : list of str, optional - Names of maps required for the metric. - - Returns - ------- - float - The magnitude of a star which has a photometric error of `crowding_error` - """ - cols = [seeing_col] units = "mag" self.crowding_error = crowding_error @@ -88,13 +92,16 @@ def __init__( super().__init__(col=cols, maps=maps, units=units, metric_name=metric_name, **kwargs) def run(self, data_slice, slice_point=None): - # Set mag_vector to the same length as starLumFunc (lower edge of mag bins) + # Set mag_vector to the same length as starLumFunc + # (lower edge of mag bins) mag_vector = slice_point[f"starMapBins_{self.filtername}"][1:] # Pull up density of stars at this point in the sky lum_func = slice_point[f"starLumFunc_{self.filtername}"] - # Calculate the crowding error using the best seeing value (in any filter?) + # Calculate the crowding error using the best seeing value + # (in any filter?) crowd_error = _comp_crowd_error(mag_vector, lum_func, seeing=min(data_slice[self.seeing_col])) - # Locate at which point crowding error is greater than user-defined limit + # Locate at which point crowding error is greater than user-defined + # limit above_crowd = np.where(crowd_error >= self.crowding_error)[0] if np.size(above_crowd) == 0: @@ -107,8 +114,29 @@ def run(self, data_slice, slice_point=None): class NstarsMetric(BaseMetric): - """Return the number of stars visible above some uncertainty limit, - taking image depth and crowding into account. + """Calculate the number of stars detectable above some uncertainty + limit, taking image depth and crowding into account. + + Parameters + ---------- + crowding_error : `float`, opt + The magnitude uncertainty from crowding in magnitudes. + Default 0.1 mags. + filtername : `str`, opt + The bandpass in which to calculate the crowding limit. Default r. + seeing_col : `str`, opt + The name of the seeing column. + m5_col : `str`, opt + The name of the m5 depth column. + maps : `list` [`str`], opt + Names of maps required for the metric. + ignore_crowding : `bool`, opt + Ignore the crowding limit. + + Returns + ------- + nstars : `float` + The number of stars above the error limit. """ def __init__( @@ -122,28 +150,6 @@ def __init__( ignore_crowding=False, **kwargs, ): - """ - Parameters - ---------- - crowding_error : float, optional - The magnitude uncertainty from crowding in magnitudes. Default 0.1 mags. - filtername: str, optional - The bandpass in which to calculate the crowding limit. Default r. - seeing_col : str, optional - The name of the seeing column. - m5_col : str, optional - The name of the m5 depth column. - maps : list of str, optional - Names of maps required for the metric. - ignore_crowding : bool (False) - Ignore the cowding limit. - - Returns - ------- - float - The number of stars above the error limit - """ - cols = [seeing_col, m5_col] units = "N stars" self.crowding_error = crowding_error @@ -157,13 +163,16 @@ def __init__( def run(self, data_slice, slice_point=None): pix_area = hp.nside2pixarea(slice_point["nside"], degrees=True) - # Set mag_vector to the same length as starLumFunc (lower edge of mag bins) + # Set mag_vector to the same length as starLumFunc + # (lower edge of mag bins) mag_vector = slice_point[f"starMapBins_{self.filtername}"][1:] # Pull up density of stars at this point in the sky lum_func = slice_point[f"starLumFunc_{self.filtername}"] - # Calculate the crowding error using the best seeing value (in any filter?) + # Calculate the crowding error using the best seeing value + # (in any filter?) crowd_error = _comp_crowd_error(mag_vector, lum_func, seeing=min(data_slice[self.seeing_col])) - # Locate at which point crowding error is greater than user-defined limit + # Locate at which point crowding error is greater than + # user-defined limit above_crowd = np.where(crowd_error >= self.crowding_error)[0] if np.size(above_crowd) == 0: @@ -171,7 +180,8 @@ def run(self, data_slice, slice_point=None): else: crowd_mag = mag_vector[max(above_crowd[0] - 1, 0)] - # Compute the coadded depth, and the mag where that depth hits the error specified + # Compute the coadded depth, and the mag where that depth + # hits the error specified coadded_depth = 1.25 * np.log10(np.sum(10.0 ** (0.8 * data_slice[self.m5_col]))) mag_limit = -2.5 * np.log10(1.0 / (self.crowding_error * (1.09 * 5))) + coadded_depth @@ -195,8 +205,17 @@ def run(self, data_slice, slice_point=None): class CrowdingMagUncertMetric(BaseMetric): - """ - Given a stellar magnitude, calculate the mean uncertainty on the magnitude from crowding. + """Calculate the mean uncertainty in magnitude due to crowding. + + Parameters + ---------- + rmag : `float` + The magnitude of the star to consider. + + Returns + ------- + mag_uncert : `float` + The uncertainty in magnitudes caused by crowding for a star of rmag. """ def __init__( @@ -209,18 +228,6 @@ def __init__( maps=["StellarDensityMap"], **kwargs, ): - """ - Parameters - ---------- - rmag : float - The magnitude of the star to consider. - - Returns - ------- - float - The uncertainty in magnitudes caused by crowding for a star of rmag. - """ - self.filtername = filtername self.seeing_col = seeing_col self.rmag = rmag diff --git a/rubin_sim/maf/metrics/cumulative_metric.py b/rubin_sim/maf/metrics/cumulative_metric.py index a18cfd97f..0d65da2e1 100644 --- a/rubin_sim/maf/metrics/cumulative_metric.py +++ b/rubin_sim/maf/metrics/cumulative_metric.py @@ -11,9 +11,9 @@ class CumulativeMetric(BaseMetric): Parameters ---------- - interp_points : `np.array` (None) - The points to interpolate the cumulative number of observations to. If None, - then the range of the data is used with a stepsize of 1. + interp_points : `np.array`, (N,) or None + The points to interpolate the cumulative number of observations to. + If None, then the range of the data is used with a stepsize of 1. """ def __init__( diff --git a/rubin_sim/maf/metrics/dcr_metric.py b/rubin_sim/maf/metrics/dcr_metric.py index 309f909fd..a77a31c0f 100644 --- a/rubin_sim/maf/metrics/dcr_metric.py +++ b/rubin_sim/maf/metrics/dcr_metric.py @@ -13,8 +13,9 @@ class DcrPrecisionMetric(BaseMetric): Parameters ---------- - atm_err : float - Minimum error in photometry centroids introduced by the atmosphere (arcseconds). Default 0.01. + atm_err : `float` + Minimum error in photometry centroids introduced by the atmosphere + (arcseconds). Default 0.01. """ def __init__( @@ -74,7 +75,8 @@ def run(self, data_slice, slice_point=None): # Now I want to compute the error if I interpolate/extrapolate to +/-1. # function is of form, y=ax. a=y/x. da = dy/x. - # Only strictly true if we know the unshifted position. But this should be a reasonable approx. + # Only strictly true if we know the unshifted position. + # But this should be a reasonable approx. slope_uncerts = position_errors / x_coord slope_uncerts2 = position_errors / x_coord2 @@ -82,10 +84,12 @@ def run(self, data_slice, slice_point=None): np.sum(1.0 / slope_uncerts**2) + np.sum(1.0 / slope_uncerts2**2) ) - # So, this will be the uncertainty in the RA or Dec offset at x= +/- 1. A.K.A., the uncertainty in the slope + # So, this will be the uncertainty in the RA or Dec offset at + # x= +/- 1. A.K.A., the uncertainty in the slope # of the line made by tan(zd)*sin(PA) vs RA offset # or the line tan(zd)*cos(PA) vs Dec offset - # Assuming we know the unshfted position of the object (or there's little covariance if we are fitting for both) + # Assuming we know the unshfted position of the object + # (or there's little covariance if we are fitting for both) result = total_slope_uncert return result diff --git a/rubin_sim/maf/metrics/exgal_m5.py b/rubin_sim/maf/metrics/exgal_m5.py index 37cd4ee51..d34af24ef 100644 --- a/rubin_sim/maf/metrics/exgal_m5.py +++ b/rubin_sim/maf/metrics/exgal_m5.py @@ -18,31 +18,40 @@ class ExgalM5(BaseMetric): Column name for five sigma depth. Default 'fiveSigmaDepth'. unit : `str`, optional Label for units. Default 'mag'. + + Returns + ------- + coadd_m5 : `float` + Coadded m5 value, corrected for galactic dust extinction. """ def __init__( self, m5_col="fiveSigmaDepth", metric_name="ExgalM5", units="mag", filter_col="filter", **kwargs ): - # Set the name for the dust map to use. This is gathered into the MetricBundle. + # Set the name for the dust map to use. + # This is gathered into the MetricBundle. maps = ["DustMap"] self.m5_col = m5_col self.filter_col = filter_col super().__init__( col=[self.m5_col, self.filter_col], maps=maps, metric_name=metric_name, units=units, **kwargs ) - # Set the default wavelength limits for the lsst filters. These are approximately correct. + # Set the default wavelength limits for the lsst filters. + # These are approximately correct. dust_properties = DustValues() self.ax1 = dust_properties.ax1 - # We will call Coaddm5Metric to calculate the coadded depth. Set it up here. + # We will call Coaddm5Metric to calculate the coadded depth. + # Set it up here. self.coaddm5_metric = Coaddm5Metric(m5_col=m5_col) def run(self, data_slice, slice_point): - """ - Compute the co-added m5 depth and then apply dust extinction to that magnitude. + """Compute the co-added m5 depth and then apply + dust extinction to that magnitude. """ m5 = self.coaddm5_metric.run(data_slice) if m5 == self.coaddm5_metric.badval: return self.badval - # Total dust extinction along this line of sight. Correct default A to this EBV value. + # Total dust extinction along this line of sight. + # Correct default A to this EBV value. a_x = self.ax1[data_slice[self.filter_col][0]] * slice_point["ebv"] return m5 - a_x diff --git a/rubin_sim/maf/metrics/galactic_plane_metrics.py b/rubin_sim/maf/metrics/galactic_plane_metrics.py index f9c6da0f9..02238f7d5 100644 --- a/rubin_sim/maf/metrics/galactic_plane_metrics.py +++ b/rubin_sim/maf/metrics/galactic_plane_metrics.py @@ -13,19 +13,24 @@ from .base_metric import BaseMetric -# These are a suite of metrics aimed at evaluating high-level quantities regarding galactic plane -# coverage. The metrics here evaluate the coverage (just number of visits and exposure time per filter) +# These are a suite of metrics aimed at evaluating high-level +# quantities regarding galactic plane coverage. +# The metrics here evaluate the coverage +# (just number of visits and exposure time per filter) # in relation to the desired coverage from the galactic plane priority map. -# There is a related metric in transientTimeSampling which evaluates the cadence weighted by this same map. +# There is a related metric in transientTimeSampling which +# evaluates the cadence weighted by this same map. TAU_OBS = np.array([2.0, 5.0, 11.0, 20.0, 46.5, 73.0]) def galplane_nvisits_thresholds(tau_obs, nyears=10): - """ "Return estimated nvisits required to well-sample lightcurves that need sampling every tau_obs (days). + """Return estimated nvisits required to well-sample lightcurves + that need sampling every tau_obs (days). - This does a very basic estimate, just counting how many visits you would have if you distributed them - at tau_obs intervals for a period of nyears, assuming a season length of 6.5 years and that visits in + This does a very basic estimate, just counting how many visits you + would have if you distributed them at tau_obs intervals for a period + of nyears, assuming a season length of 6.5 years and that visits in each night are in pairs. Parameters @@ -38,7 +43,8 @@ def galplane_nvisits_thresholds(tau_obs, nyears=10): Returns ------- n_visits_thresholds : `np.ndarray` - Estimated number of visits required to well sample lightcurves which require sampling on tau_obs + Estimated number of visits required to well sample lightcurves + which require sampling on tau_obs """ # How many nights in the survey nnights_total = 365.25 * nyears @@ -50,7 +56,8 @@ def galplane_nvisits_thresholds(tau_obs, nyears=10): def galplane_priority_map_thresholds(science_map): - """Return minimum threshold for priority maps, when considering filter balance. + """Return minimum threshold for priority maps, + when considering filter balance. Parameters ---------- @@ -84,23 +91,26 @@ def _nvisits_cut(obj, metricval): class GalPlaneFootprintMetric(BaseMetric): - """Evaluate the survey overlap with desired regions in the Galactic Plane - and Magellanic Clouds, by referencing the pre-computed priority maps provided. + """Evaluate the survey overlap with desired regions in the + Galactic Plane and Magellanic Clouds, by referencing the + pre-computed priority maps provided. These priority maps are keyed by science area (science_map) and per filter. The returned metric values are summed over all filters. Parameters ---------- science_map : `str` - Name of the priority footprint map key to use from the column headers contained in the - priority_GalPlane_footprint_map_data tables. + Name of the priority footprint map key to use from the column + headers contained in the priority_GalPlane_footprint_map_data tables. tau_obs : `np.ndarray` or `list` of `float`, opt - Timescales of minimum-required observations intervals for various classes of time variability. - Default (None), uses TAU_OBS. In general, this should be left as the default and consistent - across all galactic-plane oriented metrics. + Timescales of minimum-required observations intervals for + various classes of time variability. + Default (None), uses TAU_OBS. In general, this should be left as + the default and consistent across all galactic-plane oriented metrics. mag_cuts : `dict` of `float`, opt Magnitudes to use as cutoffs for individual image depths. - Default None uses a default set of values which correspond roughly to the 50th percentile. + Default None uses a default set of values which correspond + roughly to the 50th percentile. filter_col : `str`, opt Name of the filter column. Default 'filter'. m5_col : `str`, opt @@ -173,11 +183,14 @@ def __init__( self.reduce_order[r_name] = i + 2 def run(self, data_slice, slice_point): - """Calculate the number of observations that meet the mag_cut values at each slice_point. - Also calculate the number of observations * the priority map summed over all filter. - Return both of these values as a dictionary. + """Calculate the number of observations that meet the mag_cut values + at each slice_point. + + Also calculate the number of observations * the priority map summed + over all filter. Return both of these values as a dictionary. """ - # Check if we want to evaluate this part of the sky, or if the weight is below threshold. + # Check if we want to evaluate this part of the sky, + # or if the weight is below threshold. mapkey = gp_priority_map_components_to_keys("sum", self.science_map) priority = slice_point[mapkey] if priority <= self.priority_map_threshold: @@ -206,19 +219,23 @@ def reduce_n_obs_priority(self, metricval): class GalPlaneTimePerFilterMetric(BaseMetric): - """Evaluate the fraction of exposure time spent in each filter as a fraction of the - total exposure time dedicated to that healpix in the weighted galactic plane priority maps. + """Evaluate the fraction of exposure time spent in each filter as a + fraction of the total exposure time dedicated to that healpix in the + weighted galactic plane priority maps. Parameters ---------- scienceMap : `str` - Name of the priority footprint map key to use from the column headers contained in the + Name of the priority footprint map key to use from the column + headers contained in the priority_GalPlane_footprint_map_data tables. magCuts : `dict` of `float`, opt Magnitudes to use as cutoffs for individual image depths. - Default None uses a default set of values which correspond roughly to the 50th percentile. + Default None uses a default set of values which correspond + roughly to the 50th percentile. mjd_col : `str`, opt - Name of the observation start MJD column. Default 'observationStartMJD'. + Name of the observation start MJD column. + Default 'observationStartMJD'. exp_time_col : `str`, opt Name of the exposure time column. Default 'visitExposureTime'. filter_col : `str`, opt @@ -282,14 +299,17 @@ def __init__( def run(self, data_slice, slice_point): """Calculate the ratio of the actual on-sky exposure time per filter - compared to the ideal on-sky exposure time per filter at this point on the sky across all filters. + compared to the ideal on-sky exposure time per filter at this point + on the sky across all filters. """ - # Check if we want to evaluate this part of the sky, or if the weight is below threshold. + # Check if we want to evaluate this part of the sky, + # or if the weight is below threshold. weight_all_filters = slice_point[gp_priority_map_components_to_keys("sum", self.science_map)] if weight_all_filters <= self.priority_map_threshold: return self.badval - # Calculate the ideal weighting per filter compared to all filters at this point in the sky + # Calculate the ideal weighting per filter compared to all + # filters at this point in the sky relative_filter_weight = {} for f in self.filterlist: mapkey = gp_priority_map_components_to_keys(f, self.science_map) @@ -313,19 +333,21 @@ def run(self, data_slice, slice_point): # provided, and additional data in other filters is usually welcome exp_time_per_filter[f] = data_slice[self.exp_time_col][match].sum() - # Calculate the time on-sky in each filter that overlaps this point, and meets mag_cuts + # Calculate the time on-sky in each filter that overlaps this point, + # and meets mag_cuts total_expt_mag_cut = 0 for f in self.filterlist: total_expt_mag_cut += exp_time_per_filter[f].sum() - # normalize by the relative filter weight. Ideally metric results are close to 1. + # normalize by the relative filter weight. + # Ideally metric results are close to 1. normalized_exp_time = {} for f in self.filterlist: if total_expt_mag_cut == 0: normalized_exp_time[f] = 0 else: - # If no exposures are expected in this filter for this location, - # thais metric returns the mask val for this filter only. + # If no exposures are expected in this filter for this + # location, metric returns the mask val for this filter only. if relative_filter_weight[f] > 0: normalized_exp_time[f] = ( exp_time_per_filter[f] / total_expt_mag_cut / relative_filter_weight[f] diff --git a/rubin_sim/maf/metrics/galplane_time_sampling_metrics.py b/rubin_sim/maf/metrics/galplane_time_sampling_metrics.py index f5578efb5..f366fda50 100644 --- a/rubin_sim/maf/metrics/galplane_time_sampling_metrics.py +++ b/rubin_sim/maf/metrics/galplane_time_sampling_metrics.py @@ -1,8 +1,8 @@ -################################################################################################ +################################################################### # Metric to evaluate the transientTimeSamplingMetric # # Author - Rachel Street: rstreet@lco.global -################################################################################################ +################################################################### __all__ = ( "calc_interval_decay", "GalPlaneVisitIntervalsTimescaleMetric", @@ -48,17 +48,19 @@ class GalPlaneVisitIntervalsTimescaleMetric(BaseMetric): Parameters ---------- science_map : `str` - Name of the priority footprint map key to use from the column headers contained in the - priority_GalPlane_footprint_map_data tables. + Name of the priority footprint map key to use from the column + headers contained in the priority_GalPlane_footprint_map_data tables. tau_obs : `np.ndarray` or `list` of `float`, opt - Timescales of minimum-required observations intervals for various classes of time variability. - Default (None), uses TAU_OBS. In general, this should be left as the default and consistent - across all galactic-plane oriented metrics. + Timescales of minimum-required observations intervals for various + classes of time variability. + Default (None), uses TAU_OBS. In general, this should be left as the + default and consistent across all galactic-plane oriented metrics. mag_limit : `float`, opt Magnitude limit to use as a cutoff for various observations. Default 22.0. mjd_col : `str`, opt - The name of the observation start MJD column. Default 'observationStartMJD'. + The name of the observation start MJD column. + Default 'observationStartMJD'. m5_col : `str', opt The name of the five sigma depth column. Default 'fiveSigmaDepth'. """ @@ -80,7 +82,8 @@ def __init__( self.tau_obs = tau_obs else: self.tau_obs = TAU_OBS - # Create reduce functions for the class that are return the metric for each value in tau_obs + # Create reduce functions for the class that are return the metric + # for each value in tau_obs self.mag_limit = mag_limit self.mjd_col = mjd_col @@ -105,7 +108,8 @@ def __init__( self.reduce_order[f"reduceTau_{tau:.1f}".replace(".", "_").replace("reduce", "")] = i def run(self, data_slice, slice_point=None): - # Check if we want to evaluate this part of the sky, or if the weight is below threshold. + # Check if we want to evaluate this part of the sky, + # or if the weight is below threshold. if ( slice_point[gp_priority_map_components_to_keys("sum", self.science_map)] <= self.priority_map_threshold @@ -114,14 +118,16 @@ def run(self, data_slice, slice_point=None): # Select observations in the time sequence that fulfill the # S/N requirements: match = np.where(data_slice[self.m5_col] >= self.mag_limit)[0] - # We need at least two visits which match these requirements to calculate visit gaps + # We need at least two visits which match these requirements + # to calculate visit gaps if len(match) < 2: return self.badval # Find the time gaps between visits (in any filter) times = data_slice[self.mjd_col][match] times.sort() delta_tobs = np.diff(times) - # Compare the time gap distribution to the time gap required to characterize variability + # Compare the time gap distribution to the time gap required + # to characterize variability metric_data = {} for tau in self.tau_obs: # Normalize @@ -130,18 +136,20 @@ def run(self, data_slice, slice_point=None): class GalPlaneSeasonGapsTimescaleMetric(BaseMetric): - """Metric to evaluate the gap between sequential seasonal gaps in + """Evaluate the gap between sequential seasonal gaps in observations in a lightcurve relative to the scientifically desired sampling interval. Parameters ---------- science_map : `str` - Name of the priority footprint map key to use from the column headers contained in the + Name of the priority footprint map key to use from the column + headers contained in the priority_GalPlane_footprint_map_data tables. tau_var : `np.ndarray` or `list` of `float`, opt Timescales of variability for various classes of time variability. - Default (None), uses TAU_OBS * 5. In general, this should be left as the default and consistent + Default (None), uses TAU_OBS * 5. In general, this should be left + as the default and consistent across all galactic-plane oriented metrics. mag_limit : `float`, opt Magnitude limit to use as a cutoff for various observations. @@ -150,7 +158,8 @@ class GalPlaneSeasonGapsTimescaleMetric(BaseMetric): The typical season gap expected for a galactic plane field in days. The default, 145 days, is typical for a bulge field. mjd_col : `str`, opt - The name of the observation start MJD column. Default 'observationStartMJD'. + The name of the observation start MJD column. + Default 'observationStartMJD'. m5_col : `str', opt The name of the five sigma depth column. Default 'fiveSigmaDepth'. """ @@ -168,14 +177,16 @@ def __init__( self.science_map = science_map self.priority_map_threshold = galplane_priority_map_thresholds(self.science_map) # tau_obs is an array of minimum-required observation intervals for - # four categories of time variability; tau_var is the related timescale for the variability - # (tau_var is approximately 5 * tau_obs, in general) + # four categories of time variability; tau_var is the related timescale + # for the variability (tau_var is approximately 5*tau_obs, in general) if tau_var is not None: self.tau_var = tau_var else: self.tau_var = TAU_OBS * 5 - ### NOTE: I would recommend dropping tau_var 10 and 25 from this analysis unless the metric is changed - ### these intervals are so short they will *always* be dropped during the season gap + ### NOTE: I would recommend dropping tau_var 10 and 25 from this + # analysis unless the metric is changed + # these intervals are so short they will *always* be dropped + # during the season gap self.mag_limit = mag_limit self.expected_season_gap = expected_season_gap self.mjd_col = mjd_col @@ -194,7 +205,8 @@ def __init__( self.reduce_order[f"reduce_Tau_{tau:.1f}".replace(".", "_").replace("reduce", "")] = i def run(self, data_slice, slice_point): - # Check if we want to evaluate this part of the sky, or if the weight is below threshold. + # Check if we want to evaluate this part of the sky, + # or if the weight is below threshold. if ( slice_point[gp_priority_map_components_to_keys("sum", self.science_map)] <= self.priority_map_threshold @@ -204,17 +216,20 @@ def run(self, data_slice, slice_point): times = data_slice[self.mjd_col] times.sort() # data = np.sort(data_slice[self.mjd_col], order=self.mjd_col) - # SlicePoints ra/dec are always in radians - convert to degrees to calculate season + # SlicePoints ra/dec are always in radians - + # convert to degrees to calculate season seasons = calc_season(np.degrees(slice_point["ra"]), times) first_of_season, last_of_season = find_season_edges(seasons) - # season_lengths = times[last_of_season] - times[first_of_season] # would this match interval calc better? + # season_lengths = times[last_of_season] - times[first_of_season] + # would this match interval calc better? season_gaps = times[first_of_season][1:] - times[last_of_season][:-1] if len(season_gaps) == 0: return self.badval metric_data = {} for i, tau in enumerate(self.tau_var): metric_data[tau] = calc_interval_decay(season_gaps, tau) - # if the season gap is shorter than the expected season gap, count this as 'good' + # if the season gap is shorter than the expected season gap, + # count this as 'good' good_season_gaps = np.where(season_gaps <= self.expected_season_gap) metric_data[tau][good_season_gaps] = 1 metric_data[tau] = metric_data[tau].sum() / len(season_gaps) diff --git a/rubin_sim/maf/metrics/hourglass_metric.py b/rubin_sim/maf/metrics/hourglass_metric.py index b8aeaa57f..dd6b947ea 100644 --- a/rubin_sim/maf/metrics/hourglass_metric.py +++ b/rubin_sim/maf/metrics/hourglass_metric.py @@ -15,8 +15,10 @@ def nearest_val(A, val): class HourglassMetric(BaseMetric): - """Plot the filters used as a function of time. Must be used with the Hourglass Slicer. - Will totally fail in the arctic circle.""" + """Plot the filters used as a function of time. + Must be used with the Hourglass Slicer. + Will totally fail in the arctic circle. + """ def __init__( self, @@ -68,15 +70,18 @@ def run(self, data_slice, slice_point=None): perfilter["mjd"] = data_slice[self.mjd_col][good] perfilter["filter"] = data_slice[self.filter_col][good] - # brute force compute midnight times for all days between start and enc of data_slice + # brute force compute midnight times for all days between + # start and enc of data_slice times = Time(mjds, format="mjd") - # let's just find the midnight before and after each of the pre_night MJD values + # let's just find the midnight before and after each of the + # pre_night MJD values m_after = self.observer.midnight(times, "next") m_before = self.observer.midnight(times, "previous") midnights = np.unique(np.concatenate([m_before.mjd, m_after.mjd])) # calculating midnight can return nans? That seems bad. midnights = midnights[np.isfinite(midnights)] - # chop off any repeats. Need to round because observe.midnight values are not repeatable + # chop off any repeats. Need to round because observe.midnight + # values are not repeatable m10 = np.round(midnights * 10) _temp, indx = np.unique(m10, return_index=True) midnights = midnights[indx] @@ -105,12 +110,8 @@ def run(self, data_slice, slice_point=None): perfilter["midnight"] = midnights[indx] temp_indx = np.where(d1 < d2) perfilter["midnight"][temp_indx] = midnights[indx - 1][temp_indx] - try: - mtime = Time(pernight["midnight"], format="mjd") - except: - import pdb + mtime = Time(pernight["midnight"], format="mjd") - pdb.set_trace() pernight["twi12_rise"] = self.observer.twilight_morning_nautical(mtime, which="next").mjd pernight["twi12_set"] = self.observer.twilight_evening_nautical(mtime, which="previous").mjd diff --git a/rubin_sim/maf/metrics/long_gap_agn_metric.py b/rubin_sim/maf/metrics/long_gap_agn_metric.py deleted file mode 100644 index 17175a595..000000000 --- a/rubin_sim/maf/metrics/long_gap_agn_metric.py +++ /dev/null @@ -1,45 +0,0 @@ -__all__ = ("LongGapAGNMetric",) - -import numpy as np - -from .base_metric import BaseMetric - - -class LongGapAGNMetric(BaseMetric): - """max delta-t and average of the top-10 longest gaps.""" - - def __init__( - self, - metric_name="longGapAGNMetric", - mjdcol="observationStartMJD", - units="days", - xgaps=10, - badval=-666, - **kwargs, - ): - """Instantiate metric. - mjdcol = column name for exposure time dates - """ - cols = [mjdcol] - super(LongGapAGNMetric, self).__init__(cols, metric_name, units=units, **kwargs) - self.badval = badval - self.mjdcol = mjdcol - self.xgaps = xgaps - self.units = units - - def run(self, data_slice, slice_point=None): - metricval = np.diff(data_slice[self.mjdcol]) - return metricval - - def reduce_max_gap(self, metricval): - if metricval.size > 0: - result = np.max(metricval) - else: - result = self.badval - return result - - def reduce_average_longest_x_gaps(self, metricval): - if np.size(metricval) - self.xgaps > 0: - return np.average(np.sort(metricval)[np.size(metricval) - self.xgaps :]) - else: - return self.badval diff --git a/rubin_sim/maf/metrics/mo_metrics.py b/rubin_sim/maf/metrics/mo_metrics.py index 8d574e283..4b56c59b4 100644 --- a/rubin_sim/maf/metrics/mo_metrics.py +++ b/rubin_sim/maf/metrics/mo_metrics.py @@ -40,7 +40,8 @@ def _set_vis(sso_obs, snr_limit, snr_col, vis_col): class BaseMoMetric(BaseMetric): """Base class for the moving object metrics. - Intended to be used with the Moving Object Slicer.""" + Intended to be used with the Moving Object Slicer. + """ def __init__( self, @@ -67,7 +68,8 @@ def __init__( self.name = metric_name if self.name is None: self.name = self.__class__.__name__.replace("Metric", "", 1) - # Set badval and units, leave space for 'comment' (tied to display_dict). + # Set badval and units, leave space for 'comment' + # (tied to display_dict). self.badval = badval self.units = units self.comment = comment @@ -115,16 +117,17 @@ def run(self, sso_obs, orb, hval): Parameters ---------- - sso_obs: np.ndarray + sso_obs : `np.ndarray`, (N,) The input data to the metric (same as the parent metric). - orb: np.ndarray - The information about the orbit for which the metric is being calculated. - hval : float + orb : `np.ndarray`, (N,) + The information about the orbit for which the metric is + being calculated. + hval : `float` The H value for which the metric is being calculated. Returns ------- - float or np.ndarray or dict + metric_val : `float` or `np.ndarray` or `dict` """ raise NotImplementedError @@ -134,9 +137,10 @@ class BaseChildMetric(BaseMoMetric): Parameters ---------- - parentDiscoveryMetric: BaseMoMetric - The 'parent' metric which generated the metric data used to calculate this 'child' metric. - badval: float, optional + parentDiscoveryMetric : `~BaseMoMetric` + The 'parent' metric which generated the metric data used + calculate this 'child' metric. + badval : `float`, optional Value to return when metric cannot be calculated. """ @@ -154,18 +158,19 @@ def run(self, sso_obs, orb, hval, metric_values): Parameters ---------- - sso_obs: np.ndarray + sso_obs : `np.ndarray`, (N,) The input data to the metric (same as the parent metric). - orb: np.ndarray - The information about the orbit for which the metric is being calculated. - hval : float + orb : `np.ndarray`, (N,) + The information about the orbit for which the metric is + being calculated. + hval : `float` The H value for which the metric is being calculated. - metric_values : dict or np.ndarray + metric_values : `dict` or `np.ndarray`, (N,) The return value from the parent metric. Returns ------- - float + metric_val : `float` """ raise NotImplementedError @@ -173,13 +178,16 @@ def run(self, sso_obs, orb, hval, metric_values): class NObsMetric(BaseMoMetric): """ Count the total number of observations where an SSobject was 'visible'. + + Parameters + ---------- + snr_limit : `float` or None + If the snr_limit is None, detection of the object in a visit is + determined using the _calcVis method (completeness calculation). + If not None, the snr is calculated and used as a flat cutoff instead. """ def __init__(self, snr_limit=None, **kwargs): - """ - @ snr_limit .. if snr_limit is None, this uses the _calcVis method/completeness - if snr_limit is not None, this uses that value as a cutoff instead. - """ super().__init__(**kwargs) self.snr_limit = snr_limit @@ -193,9 +201,8 @@ def run(self, sso_obs, orb, hval): class NObsNoSinglesMetric(BaseMoMetric): - """ - Count the number of observations for an SSobject, without singles. - Don't include any observations where it was a single observation on a night. + """Count the number of observations for an SSobject, without singles. + Don't include observations where it was a single observation on a night. """ def __init__(self, snr_limit=None, **kwargs): @@ -217,10 +224,6 @@ class NNightsMetric(BaseMoMetric): """Count the number of distinct nights an SSobject is observed.""" def __init__(self, snr_limit=None, **kwargs): - """ - @ snr_limit : if SNRlimit is None, this uses _calcVis method/completeness - else if snr_limit is not None, it uses that value as a cutoff. - """ super().__init__(**kwargs) self.snr_limit = snr_limit @@ -233,7 +236,9 @@ def run(self, sso_obs, orb, hval): class ObsArcMetric(BaseMoMetric): - """Calculate the difference between the first and last observation of an SSobject.""" + """Calculate the difference in time between the first and last observation + of an SSobject. + """ def __init__(self, snr_limit=None, **kwargs): super().__init__(**kwargs) @@ -252,25 +257,29 @@ class DiscoveryMetric(BaseMoMetric): Parameters ---------- - n_obs_per_night : int, optional + n_obs_per_night : `int`, optional Number of observations required within a single night. Default 2. - t_min : float, optional + t_min : `float`, optional Minimum time span between observations in a single night, in days. Default 5 minutes (5/60/24). - t_max : float, optional + t_max : `float`, optional Maximum time span between observations in a single night, in days. Default 90 minutes. - n_nights_per_window : int, optional - Number of nights required with observations, within the track window. Default 3. - t_window : int, optional + n_nights_per_window : `int`, optional + Number of nights required with observations, within the track window. + Default 3. + t_window : `int`, optional Number of nights included in the track window. Default 15. - snr_limit : None or float, optional - SNR limit to use for observations. If snr_limit is None, (default), then it uses - the completeness calculation added to the 'vis' column (probabilistic visibility, - based on 5-sigma limit). If snr_limit is not None, it uses this SNR value as a cutoff. - metricName : str, optional + snr_limit : None or `float`, optional + SNR limit to use for observations. + If snr_limit is None, (default), then it uses + the completeness calculation added to the 'vis' column + (probabilistic visibility, based on 5-sigma limit). + If snr_limit is not None, it uses this SNR value as a cutoff. + metricName : `str`, optional The metric name to use. - Default will be to construct Discovery_nObsPerNightxnNightsPerWindowintWindow. + Default will be to construct + Discovery_nObsPerNightxnNightsPerWindowintWindow. """ def __init__( @@ -356,8 +365,9 @@ def run(self, sso_obs, orb, hval): tidx = np.where((dtimes >= self.t_min) & (dtimes <= self.t_max))[0] if len(tidx) > 0: good[c] = 1 - # 'good' provides mask for observations which could count as 'good to make tracklets' - # against sso_obs[vis_sort][n_idx_many]. Now identify tracklets which can make tracks. + # 'good' provides mask for observations which could count as + # 'good to make tracklets' against sso_obs[vis_sort][n_idx_many]. + # Now identify tracklets which can make tracks. good_idx = vis_sort[n_idx_many][good == 1] good_idx_ends = vis_sort[n_idx_many_end][good == 1] # print 'good tracklets', nights[good_idx] @@ -367,7 +377,8 @@ def run(self, sso_obs, orb, hval): np.roll(sso_obs[self.night_col][vis][good_idx], 1 - self.n_nights_per_window) - sso_obs[self.night_col][vis][good_idx] ) - # Identify the index in sso_obs[vis][good_idx] (sorted by mjd) where the discovery opportunity starts. + # Identify the index in sso_obs[vis][good_idx] (sorted by mjd) + # where the discovery opportunity starts. start_idxs = np.where((delta_nights >= 0) & (delta_nights <= self.t_window))[0] # Identify the index where the discovery opportunity ends. end_idxs = np.zeros(len(start_idxs), dtype="int") @@ -391,15 +402,13 @@ def run(self, sso_obs, orb, hval): class DiscoveryNChancesMetric(BaseChildMetric): """Calculate total number of discovery opportunities for an SSobject. - Retirms total number of discovery opportunities. + Returns total number of discovery opportunities. Child metric to be used with the Discovery Metric. """ def __init__( self, parent_discovery_metric, - # night_start=None, - # night_end=None, badval=0, **kwargs, ): @@ -407,7 +416,8 @@ def __init__( self.night_start = None # night_start self.night_end = None # night_end self.snr_limit = parent_discovery_metric.snr_limit - # Update the metric name to use the night_start/night_end values, if an overriding name is not given. + # Update the metric name to use the night_start/night_end values, + # if an overriding name is not given. if "metric_name" not in kwargs: if self.night_start is not None: self.name = self.name + "_n%d" % (self.night_start) @@ -415,34 +425,16 @@ def __init__( self.name = self.name + "_n%d" % (self.night_end) def run(self, sso_obs, orb, hval, metric_values): - """Return the number of different discovery chances we had for each object/H combination.""" - return metric_values["n_chances"] - """ - vis = _set_vis(sso_obs, self.snr_limit, self.snr_col, self.vis_col) - if len(vis) == 0: - return self.badval - if self.night_start is None and self.night_end is None: - return len(metric_values["start"]) - # Otherwise, we have to sort out what night the discovery chances happened on. - vis_sort = np.argsort(sso_obs[self.mjd_col][vis]) - nights = sso_obs[self.night_col][vis][vis_sort] - start_nights = nights[metric_values["start"]] - end_nights = nights[metric_values["end"]] - if self.night_end is None and self.night_start is not None: - valid = np.where(start_nights >= self.night_start)[0] - elif self.night_start is None and self.night_end is not None: - valid = np.where(end_nights <= self.night_end)[0] - else: - # And we only end up here if both were not None. - valid = np.where( - (start_nights >= self.night_start) & (end_nights <= self.night_end) - )[0] - return len(valid) + """Return the number of different discovery chances we + had for each object/H combination. """ + return metric_values["n_chances"] class DiscoveryNObsMetric(BaseChildMetric): - """Calculates the number of observations in the first discovery track of an SSobject.""" + """Calculates the number of observations in the first discovery + track of an SSobject. + """ def __init__(self, parent_discovery_metric, badval=0, **kwargs): super().__init__(parent_discovery_metric, badval=badval, **kwargs) @@ -527,7 +519,8 @@ def run(self, sso_obs, orb, hval, metric_values): class DiscoveryEclonlatMetric(BaseChildMetric): - """Returns the ecliptic lon/lat and solar elong of the first discovery track of an SSobject.""" + """Returns the ecliptic lon/lat and solar elong of the first discovery + track of an SSobject.""" def __init__(self, parent_discovery_metric, badval=None, **kwargs): super().__init__(parent_discovery_metric, badval=badval, **kwargs) @@ -574,7 +567,8 @@ class ActivityOverTimeMetric(BaseMoMetric): Counts the time periods where we would have a chance to detect activity on a moving object. - Splits observations into time periods set by 'window', then looks for observations within each window, + Splits observations into time periods set by 'window', + then looks for observations within each window, and reports what fraction of the total windows receive 'nObs' visits. """ @@ -590,7 +584,8 @@ def __init__(self, window, snr_limit=5, survey_years=10.0, metric_name=None, **k self.units = "%.1f Day Windows" % (self.window) def run(self, sso_obs, orb, hval): - # For cometary activity, expect activity at the same point in its orbit at the same time, mostly + # For cometary activity, expect activity at the same point in its + # orbit at the same time, mostly # For collisions, expect activity at random times vis = _set_vis(sso_obs, self.snr_limit, self.snr_col, self.vis_col) if len(vis) == 0: @@ -601,7 +596,8 @@ def run(self, sso_obs, orb, hval): class ActivityOverPeriodMetric(BaseMoMetric): - """Count fraction of object period we could identify activity for an SSobject. + """Count fraction of object period we could identify activity + for an SSobject. Count the fraction of the orbit (when split into n_bins) that receive observations, in order to have a chance to detect activity. @@ -638,7 +634,8 @@ def __init__( self.units = "%.1f deg" % (np.degrees(self.bin_size)) def run(self, sso_obs, orb, hval): - # For cometary activity, expect activity at the same point in its orbit at the same time, mostly + # For cometary activity, expect activity at the same point in its + # orbit at the same time, mostly # For collisions, expect activity at random times if self.a_col in orb.keys(): a = orb[self.a_col] @@ -667,15 +664,21 @@ def run(self, sso_obs, orb, hval): class MagicDiscoveryMetric(BaseMoMetric): - """Count the number of nights with discovery opportunities with very good software for an SSobject.""" + """Count the number of nights with discovery opportunities + with very good software for an SSobject. + + Parameters + ---------- + n_obs : `int`, opt + Total number of observations required for discovery. + t_window : `float`, opt + The timespan of the discovery window (days). + snr_limit : `float` or None + If None, uses the probabilistic detection likelihood. + If float, uses the SNR value as a flat cutoff value. + """ def __init__(self, n_obs=6, t_window=60, snr_limit=None, **kwargs): - """ - @ n_obs = the total number of observations required for 'discovery' - @ t_window = the timespan of the discovery window. - @ snr_limit .. if snr_limit is None then uses 'completeness' calculation, - .. if snr_limit is not None, then uses this value as a cutoff. - """ super().__init__(**kwargs) self.snr_limit = snr_limit self.n_obs = n_obs @@ -683,7 +686,6 @@ def __init__(self, n_obs=6, t_window=60, snr_limit=None, **kwargs): self.badval = 0 def run(self, sso_obs, orb, hval): - """SsoObs = Dataframe, orb=Dataframe, hval=single number.""" # Calculate visibility for this orbit at this H. vis = _set_vis(sso_obs, self.snr_limit, self.snr_col, self.vis_col) if len(vis) < self.n_obs: @@ -699,16 +701,22 @@ def run(self, sso_obs, orb, hval): class HighVelocityMetric(BaseMoMetric): """Count number of times an SSobject appears trailed. - Count the number of times an asteroid is observed with a velocity high enough to make it appear - trailed by a factor of (psf_factor)*PSF - i.e. velocity >= psf_factor * seeing / visitExpTime. + Count the number of times an asteroid is observed with a velocity + high enough to make it appear trailed by a factor of (psf_factor)*PSF - + i.e. velocity >= psf_factor * seeing / visitExpTime. Simply counts the total number of observations with high velocity. + + Parameters + ---------- + psf_factor : `float`, opt + The factor to multiply the seeing/VisitExpTime by to compare against + velocity. + snr_limit : `float` or None + If None, uses the probabilistic detection likelihood. + If float, uses the SNR value as a flat cutoff value. """ def __init__(self, psf_factor=2.0, snr_limit=None, velocity_col="velocity", **kwargs): - """ - @ psf_factor = factor to multiply seeing/visitExpTime by - (velocity(deg/day) >= 24*psf_factor*seeing(")/visitExptime(s)) - """ super().__init__(**kwargs) self.velocity_col = velocity_col self.snr_limit = snr_limit @@ -727,29 +735,35 @@ def run(self, sso_obs, orb, hval): class HighVelocityNightsMetric(BaseMoMetric): - """Count the number of discovery opportunities (via trailing) for an SSobject. + """Count the number of discovery opportunities (via trailing) for an + SSobject. - Determine the first time an asteroid is observed is observed with a velocity high enough to make - it appear trailed by a factor of psf_factor*PSF with n_obs_per_night observations within a given night. + Determine the first time an asteroid is observed is observed with a + velocity high enough to make it appear trailed by a factor of + psf_factor*PSF with n_obs_per_night observations within a given night. Parameters ---------- - psf_factor: float, optional - Object velocity (deg/day) must be >= 24 * psf_factor * seeingGeom (") / visitExpTime (s). + psf_facto r: `float`, optional + Object velocity (deg/day) must be + >= 24 * psf_factor * seeingGeom (") / visitExpTime (s). Default is 2 (i.e. object trailed over 2 psf's). - n_obs_per_night: int, optional + n_obs_per_night : `int`, optional Number of observations per night required. Default 2. - snr_limit: float or None - If snr_limit is set as a float, then requires object to be above snr_limit SNR in the image. - If snr_limit is None, this uses the probabilistic 'visibility' calculated by the vis stacker, - which means SNR ~ 5. Default is None. - velocity_col: str, optional - Name of the velocity column in the obs file. Default 'velocity'. (note this is deg/day). + snr_limit : `float` or None + If snr_limit is set as a float, then requires object to be above + snr_limit SNR in the image. + If snr_limit is None, this uses the probabilistic 'visibility' + calculated by the vis stacker, which means SNR ~ 5. + Default is None. + velocity_col : `str`, optional + Name of the velocity column in the obs file. + Default 'velocity'. (note this is deg/day). Returns ------- - float - The time of the first detection where the conditions are satisifed. + time : `float` + The time of the first detection where the conditions are satisfed. """ def __init__(self, psf_factor=2.0, n_obs_per_night=2, snr_limit=None, velocity_col="velocity", **kwargs): @@ -788,45 +802,54 @@ def run(self, sso_obs, orb, hval): class LightcurveInversionAsteroidMetric(BaseMoMetric): - """ - This metric is generally applicable to NEOs and MBAs - inner solar system objects. + """Evaluate the liklihood that the detections could be used to enable + lightcurve inversion. This metric is generally applicable only to inner + solar system objects (NEOs, MBAs). - Determine if the cumulative sum of observations of a target are enough to enable lightcurve - inversion for shape modeling. For this to be true, multiple conditions need to be + Determine if the cumulative sum of observations of a target are + enough to enable lightcurve inversion for shape modeling. + For this to be true, multiple conditions need to be satisfied: - 1) The SNR-weighted number of observations (each observation is weighted by its SNR, up to a max of 100) - must be larger than the threshhold weight_det (default 50) - 2) Ecliptic longitudinal coverage needs to be at least 90 degrees, and the absolute deviation - needs to be at least 1/8th the longitudinal coverage. + 1) The SNR-weighted number of observations (each observation is weighted + by its SNR, up to a max of 100) must be larger than the + threshold weight_det (default 50) + 2) Ecliptic longitudinal coverage needs to be at least 90 degrees, + and the absolute deviation needs to be at least 1/8th the + longitudinal coverage. 3) The phase angle coverage needs to span at least 5 degrees. - For evaluation of condition 2, the median ecliptic longitude is subtracted from all longitudes, - and the modulo 360 of those values is taken. This ensures that the wrap around 360 is handled - correctly. + For evaluation of condition 2, the median ecliptic longitude is + subtracted from all longitudes, and the modulo 360 of those values + is taken. This ensures that the wrap around 360 is handled correctly. For more information on the above conditions, please see https://docs.google.com/document/d/1GAriM7trpTS08uanjUF7PyKALB2JBTjVT7Y6R30i0-8/edit?usp=sharing - Contributed by Steve Chesley, Wes Fraser, Josef Durech, and the inner solar system working group. + + Contributed by Steve Chesley, Wes Fraser, Josef Durech, and the + inner solar system working group. Parameters ---------- - weight_det: float, optional - The SNR-weighted number of detections required (per bandpass in any ONE of the filters in filterlist). + weight_det : `float`, optional + The SNR-weighted number of detections required (per bandpass in any + ONE of the filters in filterlist). Default 50. - snr_limit: float or None, optional - If snr_limit is set as a float, then requires object to be above snr_limit SNR in the image. - If snr_limit is None, this uses the probabilistic 'visibility' calculated by the vis stacker, + snr_limit : `float` or None, optional + If snr_limit is set as a float, then requires object to be + above snr_limit SNR in the image. + If snr_limit is None, this uses the probabilistic 'visibility' + calculated by the vis stacker, which means SNR ~ 5. Default is None. - snr_max: float, optional + snr_max : `float`, optional Maximum value toward the SNR-weighting to consider. Default 100. - filterlist: list of str, optional - The filters which the lightcurve inversion could be based on. Requirements must be met in one of - these filters. + filterlist : `list` [`str`], optional + The filters which the lightcurve inversion could be based on. + Requirements must be met in one of these filters. Returns ------- - int + metric_value : `int` 0 (could not perform lightcurve inversion) or 1 (could) """ @@ -856,9 +879,11 @@ def run(self, sso_obs, orb, hval): match = np.where(sso_obs[self.filter_col] == f) snr_sum = np.sum(clip_snr[match]) / self.snr_max if snr_sum < self.weight_det: - # Do not have enough SNR-weighted observations, so skip on to the next filter. + # Do not have enough SNR-weighted observations, + # so skip on to the next filter. continue - # Is the ecliptic longitude coverage for the visible observations sufficient? + # Is the ecliptic longitude coverage for the visible + # observations sufficient? # Is the phase coverage sufficient? vis = np.where(clip_snr[match] > 0) ec_l = sso_obs["ecLon"][match][vis] @@ -869,7 +894,8 @@ def run(self, sso_obs, orb, hval): d_l = np.max(ec_l) - np.min(ec_l) # Calculate the range of the phase angle dp = np.max(phase_angle) - np.min(phase_angle) - # Metric requirement is that d_l >= 90 deg, absolute deviation is greater than d_l/8 + # Metric requirement is that d_l >= 90 deg, absolute + # deviation is greater than d_l/8 # and then that the phase coverage is more than 5 degrees. # Stop as soon as find a case where this is true. if d_l >= 90.0 and a_dev >= d_l / 8 and dp >= 5: @@ -879,40 +905,50 @@ def run(self, sso_obs, orb, hval): class ColorAsteroidMetric(BaseMoMetric): - """ - This metric is appropriate for MBAs and NEOs, and other inner solar system objects. - - The metric evaluates if the SNR-weighted number of observations are enough to - determine an approximate lightcurve and phase function -- and from this, - then a color for the asteroid can be determined. - The assumption is that you must fit the lightcurve/phase function in each bandpass, - and could do this well-enough if you have at least weight_det SNR-weighted observations - in the bandpass. - e.g. to find a g-r color, you must have 10 (SNR-weighted) obs in g and 10 in r. + """Calculate the likelihood of being able to calculate the color of an + object. This metric is appropriate for MBAs and NEOs, + and other inner solar system objects. + + The metric evaluates if the SNR-weighted number of observations are + enough to determine an approximate lightcurve and phase function -- + and from this, then a color for the asteroid can be determined. + The assumption is that you must fit the lightcurve/phase function + in each bandpass, and could do this well-enough if you have at least + weight_det SNR-weighted observations in the bandpass. + e.g. to find a g-r color, you must have 10 (SNR-weighted) obs in g + and 10 in r. For more details, see https://docs.google.com/document/d/1GAriM7trpTS08uanjUF7PyKALB2JBTjVT7Y6R30i0-8/edit?usp=sharing - Contributed by Wes Fraser, Steven Chesley & the inner solar system working group. + + Contributed by Wes Fraser, Steven Chesley + & the inner solar system working group. Parameters ---------- weight_det: float, optional - The SNR-weighted number of detections required (per bandpass in any ONE of the filters in filterlist). + The SNR-weighted number of detections required (per bandpass in any + ONE of the filters in filterlist). Default 10. snr_limit: float or None, optional - If snr_limit is set as a float, then requires object to be above snr_limit SNR in the image. - If snr_limit is None, this uses the probabilistic 'visibility' calculated by the vis stacker, + If snr_limit is set as a float, then requires object to be above + snr_limit SNR in the image. + If snr_limit is None, this uses the probabilistic 'visibility' + calculated by the vis stacker, which means SNR ~ 5. Default is None. snr_max: float, optional Maximum value toward the SNR-weighting to consider. Default 20. Returns ------- - int - An integer 'flag' that indicates whether the mean magnitude (and thus a color) was determined in: + flag : `int` + An integer 'flag' that indicates whether the mean magnitude + (and thus a color) was determined in: 0 = no bands - 1 = g and (r or i) and (z or y). i.e. obtain colors g-r or g-i PLUS g-z or g-y - 2 = Any 4 different filters (from grizy). i.e. colors = g-r, r-i, i-z, OR r-i, i-z, z-y.. + 1 = g and (r or i) and (z or y). + i.e. obtain colors g-r or g-i PLUS g-z or g-y + 2 = Any 4 different filters (from grizy). + i.e. colors = g-r, r-i, i-z, OR r-i, i-z, z-y.. 3 = All 5 from grizy. i.e. colors g-r, r-i, i-z, z-y. 4 = All 6 filters (ugrizy) -- best possible! add u-g. """ @@ -944,14 +980,17 @@ def run(self, sso_obs, orb, hval): # Now assign a flag: # 0 = no bands - # 1 = g and (r or i) and (z or y). i.e. obtain colors g-r or g-i PLUS g-z or g-y - # 2 = Any 4 different filters (from grizy). i.e. colors = g-r, r-i, i-z, OR r-i, i-z, z-y.. + # 1 = g and (r or i) and (z or y). + # i.e. obtain colors g-r or g-i PLUS g-z or g-y + # 2 = Any 4 different filters (from grizy). + # i.e. colors = g-r, r-i, i-z, OR r-i, i-z, z-y.. # 3 = All 5 from grizy. i.e. colors g-r, r-i, i-z, z-y. # 4 = All 6 filters (ugrizy) -- best possible! add u-g. all_six = set(self.filterlist) good_five = set(["g", "r", "i", "z", "y"]) - if len(filter_weight) == 0: # this lets us stop evaluating here if possible. + if len(filter_weight) == 0: + # this lets us stop evaluating here if possible. flag = 0 elif all_six.intersection(filter_weight) == all_six: flag = 4 @@ -974,38 +1013,50 @@ def run(self, sso_obs, orb, hval): class LightcurveColorOuterMetric(BaseMoMetric): - """ - This metric is appropriate for outer solar system objects, such as TNOs and SDOs. + """Calculate the liklihood of being able to calculate a color and + lightcurve for outer solar system objects. - This metric evaluates whether the number of observations is sufficient to fit a lightcurve - in a primary and secondary bandpass. The primary bandpass requires more observations than - the secondary. Essentially, it's a complete lightcurve in one or both bandpasses, with at + This metric is appropriate for outer solar system objects, + such as TNOs and SDOs. + + This metric evaluates whether the number of observations is + sufficient to fit a lightcurve in a primary and secondary bandpass. + The primary bandpass requires more observations than the secondary. + Essentially, it's a complete lightcurve in one or both bandpasses, with at least a semi-complete lightcurve in the secondary band. - The lightcurve/color can be calculated with any two of the bandpasses in filterlist. + The lightcurve/color can be calculated with any two of the + bandpasses in filterlist. + Contributed by Wes Fraser. Parameters ---------- - snr_limit: float or None, optional - If snr_limit is set as a float, then requires object to be above snr_limit SNR in the image. - If snr_limit is None, this uses the probabilistic 'visibility' calculated by the vis stacker, + snr_limit : `float` or None, optional + If snr_limit is set as a float, then requires object to be above + snr_limit SNR in the image. + If snr_limit is None, this uses the probabilistic 'visibility' + calculated by the vis stacker, which means SNR ~ 5. Default is None. - num_req: int, optional + num_req : `int`, optional Number of observations required for a lightcurve fitting. Default 30. - num_sec_filt: int, optional - Number of observations required in a secondary band for color only. Default 20. - filterlist: list of str, optional + num_sec_filt : `int`, optional + Number of observations required in a secondary band for color only. + Default 20. + filterlist : `list` [`str`], optional Filters that the primary/secondary measurements can be in. Returns ------- - int + flag : `ont` A flag that indicates whether a color/lightcurve was generated in: - 0 = no lightcurve (although may have had 'color' in one or more band) - 1 = a lightcurve in a single filter (but no additional color information) + 0 = no lightcurve + (although may have had 'color' in one or more band) + 1 = a lightcurve in a single filter + (but no additional color information) 2+ = lightcurves in more than one filter (or lightcurve + color) - e.g. lightcurve in 2 bands, with additional color information in another = 3. + e.g. lightcurve in 2 bands, + with additional color information in another = 3. """ def __init__( @@ -1046,31 +1097,35 @@ def run(self, sso_obs, orb, hval): class InstantaneousColorMetric(BaseMoMetric): - """Identify SSobjects which could have observations suitable to determine colors. + """Identify SSobjects which could have observations suitable to + determine instanteous colors. - Generally, this is not the mode LSST would work in - the lightcurves of the objects - mean that the time interval would have to be quite short. + Generally, this is not the mode LSST would work in - + the lightcurves of the objects mean that the time interval would have to + be quite short. - This is roughly defined as objects which have more than n_pairs pairs of observations - with SNR greater than snr_limit, in bands bandOne and bandTwo, within n_hours. + This is roughly defined as objects which have more than n_pairs pairs + of observations with SNR greater than snr_limit, + in bands bandOne and bandTwo, within n_hours. Parameters ---------- - n_pairs: int, optional - The number of pairs of observations (in each band) that must be within n_hours - Default 1 - snr_limit: float, optional + n_pairs : `int`, optional + The number of pairs of observations (in each band) that must be + within n_hours. Default 1. + snr_limit : `float`, optional The SNR limit for the observations. Default 10. - n_hours: float, optional - The time interval between observations in the two bandpasses (hours). Default 0.5 hours. - b_one: str, optional + n_hours : `float`, optional + The time interval between observations in the two bandpasses (hours). + Default 0.5 hours. + b_one : `str`, optional The first bandpass for the color. Default 'g'. - b_two: str, optional + b_two : `str`, optional The second bandpass for the color. Default 'r'. Returns ------- - int + flag : `int` 0 (no color possible under these constraints) or 1 (color possible). """ @@ -1117,50 +1172,64 @@ def run(self, sso_obs, orb, hval): class KnownObjectsMetric(BaseMoMetric): - """Identify SSobjects which could be classified as 'previously known' based on their peak V magnitude. - This is most appropriate for NEO surveys, where most of the sky has been covered so the exact location + """Identify SSobjects which could be classified as 'previously known' + based on their peak V magnitude. + This is most appropriate for NEO surveys, where most of the sky has + been covered so the exact location (beyond being in the visible sky) is not as important. Default parameters tuned to match NEO survey capabilities. Returns the time at which each first reached that threshold V magnitude. - The default values are calibrated using the NEOs larger than 140m discovered in the last 20 years - and assuming a 30% completeness in 2017. + The default values are calibrated using the NEOs larger than 140m + discovered in the last 20 years and assuming a 30% completeness in 2017. + + Note: the default parameteres here were set up in ~2012, and are likely + out of date (potentially adding another epoch of discovery). Parameters ----------- - elong_thresh : float, optional - The cutoff in solar elongation to consider an object 'visible'. Default 100 deg. - v_mag_thresh1 : float, optional + elong_thresh : `float`, optional + The cutoff in solar elongation to consider an object 'visible'. + Default 100 deg. + v_mag_thresh1 : `float`, optional The magnitude threshold for previously known objects. Default 20.0. - eff1 : float, optional + eff1 : `float`, optional The likelihood of actually achieving each individual input observation. - If the input observations include one observation per day, an 'eff' value of 0.3 would - mean that (on average) only one third of these observations would be achieved. - This is similar to the level for LSST, which can cover the visible sky every 3-4 days. + If the input observations include one observation per day, + an 'eff' value of 0.3 would mean that (on average) only one third + of these observations would be achieved. This is similar to the level + for LSST, which can cover the visible sky every 3-4 days. Default 0.1 - t_switch1 : float, optional - The (MJD) time to switch between v_mag_thresh1 + eff1 to v_mag_thresh2 + eff2, e.g. - the end of the first period. + t_switch1 : `float`, optional + The (MJD) time to switch between v_mag_thresh1 + eff1 to + v_mag_thresh2 + eff2, e.g. the end of the first period. Default 53371 (2005). - v_mag_thresh2 : float, optional + v_mag_thresh2 : `float`, optional The magnitude threshhold for previously known objects. Default 22.0. - This is based on assuming PS and other surveys will be efficient down to V=22. - eff2 : float, optional - The efficiency of observations during the second period of time. Default 0.1 - t_switch2 : float, optional - The (MJD) time to switch between v_mag_thresh2 + eff2 to v_mag_thresh3 + eff3. + This is based on assuming PS and other surveys will be efficient + down to V=22. + eff2 : `float`, optional + The efficiency of observations during the second period of time. + Default 0.1 + t_switch2 : `float`, optional + The (MJD) time to switch between v_mag_thresh2 + eff2 to + v_mag_thresh3 + eff3. Default 57023 (2015). - v_mag_thresh3 : float, optional - The magnitude threshold during the third period. Default 22.0, based on PS1 + Catalina. - eff3 : float, optional + v_mag_thresh3 : `float`, optional + The magnitude threshold during the third period. + Default 22.0, based on PS1 + Catalina. + eff3 : `float`, optional The efficiency of observations during the third period. Default 0.1 - t_switch3 : float, optional - The (MJD) time to switch between v_mag_thresh3 + eff3 to v_mag_thresh4 + eff4. + t_switch3 : `float`, optional + The (MJD) time to switch between v_mag_thresh3 + eff3 + to v_mag_thresh4 + eff4. Default 59580 (2022). - v_mag_thresh4 : float, optional - The magnitude threshhold during the fourth (last) period. Default 22.0, based on PS1 + Catalina. - eff4 : float, optional - The efficiency of observations during the fourth (last) period. Default 0.2 + v_mag_thresh4 : `float`, optional + The magnitude threshhold during the fourth (last) period. + Default 22.0, based on PS1 + Catalina. + eff4 : `float`, optional + The efficiency of observations during the fourth (last) period. + Default 0.2 """ def __init__( @@ -1168,7 +1237,7 @@ def __init__( elong_thresh=100.0, v_mag_thresh1=20.0, eff1=0.1, - t_switch1=53371, # XXX--maybe swap to survey_start_mjd and then delta_t's + t_switch1=53371, v_mag_thresh2=21.5, eff2=0.1, t_switch2=57023, diff --git a/rubin_sim/maf/metrics/pair_metric.py b/rubin_sim/maf/metrics/pair_metric.py index d282699a4..032f1180c 100644 --- a/rubin_sim/maf/metrics/pair_metric.py +++ b/rubin_sim/maf/metrics/pair_metric.py @@ -6,8 +6,17 @@ class PairMetric(BaseMetric): - """ - Count the number of pairs that could be used for Solar System object detection + """Count the number of pairs of visits that could be used for + Solar System object detection. + + Parameters + ---------- + match_min : `float` + Minutes after first observation to count something as a match. + match_max : `float` + Minutes after first observation to count something as a match. + bin_size : `float` + bin_size to use (minutes). """ def __init__( @@ -19,16 +28,6 @@ def __init__( bin_size=5.0, **kwargs, ): - """ - Parameters - ---------- - match_min : float (20.) - Minutes after first observation to count something as a match - match_max : float (40.) - Minutes after first observation to count something as a match - bin_size : float (5.) - bin_size to use (minutes) - """ self.mjd_col = mjd_col self.bin_size = bin_size / 60.0 / 24.0 self.match_min = match_min / 60.0 / 24.0 @@ -47,7 +46,8 @@ def run(self, data_slice, slice_point=None): nbin_max = np.round(self.match_max / self.bin_size) bins_to_check = np.arange(nbin_min, nbin_max + 1, 1) bins_w_obs = np.where(hist > 0)[0] - # now, for each bin with an observation, need to check if there is a bin + # now, for each bin with an observation, + # need to check if there is a bin # far enough ahead that is also populated. result = 0 for binadd in bins_to_check: diff --git a/rubin_sim/maf/metrics/periodic_detect_metric.py b/rubin_sim/maf/metrics/periodic_detect_metric.py index 445bfe5e7..7cb4cbd42 100644 --- a/rubin_sim/maf/metrics/periodic_detect_metric.py +++ b/rubin_sim/maf/metrics/periodic_detect_metric.py @@ -9,30 +9,38 @@ class PeriodicDetectMetric(BaseMetric): - """Determine if we would be able to classify an object as periodic/non-uniform, using an F-test - The idea here is that if a periodic source is aliased, it will be indistinguishable from a constant source, - so we can find a best-fit constant, and if the reduced chi-squared is ~1, we know we are aliased. + """Determine if we would be able to classify an object as + periodic/non-uniform, using an F-test. + + The idea here is that if a periodic source is aliased, it will be + indistinguishable from a constant source, + so we can find a best-fit constant, and if the reduced chi-squared is ~1, + we know we are aliased. Parameters ---------- - - period : float (2) or array - The period of the star (days). Can be a single value, or an array. If an array, amplitude and starMag + period : `float` or `array` + The period of the star (days). + Can be a single value, or an array. If an array, amplitude and starMag should be arrays of equal length. - amplitude : floar (0.1) - The amplitude of the stellar variablility (mags). - starMag : float (20.) + amplitude : `float` + The amplitude of the stellar variability, (mags). + starMag : `float` The mean magnitude of the star in r (mags). - sig_level : float (0.05) - The value to use to compare to the p-value when deciding if we can reject the null hypothesis. - sed_template : str ('F') - The stellar SED template to use to generate realistic colors (default is an F star, so RR Lyrae-like) + sig_level : `float` + The value to use to compare to the p-value when deciding + if we can reject the null hypothesis. + sed_template : `str` + The stellar SED template to use to generate realistic colors + (default is an F star, so RR Lyrae-like) Returns ------- - - 1 if we would detect star is variable, 0 if it is well-fit by a constant value. If using arrays to test multiple - period-amplitude-mag combinations, will be the sum of the number of detected stars. + flag : `int` + Returns 1 if we would detect star is variable, + 0 if it is well-fit by a constant value. + If using arrays to test multiple period-amplitude-mag combinations, + will be the sum of the number of detected stars. """ def __init__( @@ -53,7 +61,8 @@ def __init__( self.filter_col = filter_col if np.size(periods) == 1: self.periods = [periods] - # Using the same magnitude for all filters. Could expand to fit the mean in each filter. + # Using the same magnitude for all filters. + # Could expand to fit the mean in each filter. self.star_mags = [star_mags] self.amplitudes = [amplitudes] else: @@ -99,11 +108,14 @@ def run(self, data_slice, slice_point=None): weights = 1.0 / (delta_m**2) weighted_mean = np.sum(weights * lc) / np.sum(weights) chi_sq_1 += np.sum(((lc - weighted_mean) ** 2 / delta_m**2)) - # Yes, I'm fitting magnitudes rather than flux. At least I feel kinda bad about it. - # F-test for nested models Regression problems: https://en.wikipedia.org/wiki/F-test + # Yes, I'm fitting magnitudes rather than flux. + # At least I feel kinda bad about it. + # F-test for nested models Regression problems: + # https://en.wikipedia.org/wiki/F-test f_numerator = (chi_sq_1 - chi_sq_2) / (p2 - p1) f_denom = 1.0 - # This is just reduced chi-squared for the more complicated model, so should be 1. + # This is just reduced chi-squared for the more + # complicated model, so should be 1. f_val = f_numerator / f_denom # Has DoF (p2-p1, n-p2) # https://stackoverflow.com/questions/21494141/how-do-i-do-a-f-test-in-python/21503346 diff --git a/rubin_sim/maf/metrics/phase_gap_metric.py b/rubin_sim/maf/metrics/phase_gap_metric.py index 8ae527189..8ec0142f1 100644 --- a/rubin_sim/maf/metrics/phase_gap_metric.py +++ b/rubin_sim/maf/metrics/phase_gap_metric.py @@ -8,21 +8,26 @@ class PhaseGapMetric(BaseMetric): - """ - Measure the maximum gap in phase coverage for observations of periodic variables. + """Measure the maximum gap in phase coverage for + observations of periodic variables. Parameters ---------- - col: str, optional + col : `str`, optional Name of the column to use for the observation times (MJD) - n_periods: int, optional + n_periods : `int`, optional Number of periods to test - period_min: float, optional + period_min : `float`, optional Minimum period to test, in days. - period_max: float, optional + period_max : `float`, optional Maximum period to test, in days - n_visits_min: int, optional + n_visits_min : `int`, optional Minimum number of visits necessary before looking for the phase gap. + + Returns + ------- + metric_value : `dict` {`periods`: `float`, `maxGaps` : `float`} + Calculates a dictionary of max gap in phase coverage for each period. """ def __init__( @@ -91,8 +96,11 @@ def reduce_largest_gap(self, metric_val): return np.max(metric_val["maxGaps"]) -# To fit a periodic source well, you need to cover the full phase, and fit the amplitude. +# To fit a periodic source well, you need to cover the full phase, +# and fit the amplitude. class PeriodicQualityMetric(BaseMetric): + """Evaluate phase coverage over a given period. + """ def __init__( self, mjd_col="observationStartMJD", @@ -111,7 +119,8 @@ def __init__( ) def _calc_phase(self, data_slice): - """1 is perfectly balanced phase coverage, 0 is no effective coverage.""" + """1 is perfectly balanced phase coverage, 0 is no effective coverage. + """ angles = data_slice[self.mjd_col] % self.period angles = angles / self.period * 2.0 * np.pi x = np.cos(angles) @@ -125,7 +134,9 @@ def _calc_phase(self, data_slice): return 1.0 - vector_off def _calc_amp(self, data_slice): - """Fractional SNR on the amplitude, testing for a variety of possible phases""" + """Fractional SNR on the amplitude, + testing for a variety of possible phases. + """ phases = np.arange(0, np.pi, np.pi / 8.0) snr = m52snr(self.star_mag, data_slice[self.m5_col]) amp_snrs = np.sin(data_slice[self.mjd_col] / self.period * 2 * np.pi + phases[:, np.newaxis]) * snr diff --git a/rubin_sim/maf/metrics/qso_number_counts_metric.py b/rubin_sim/maf/metrics/qso_number_counts_metric.py index be056a5da..6900c95aa 100644 --- a/rubin_sim/maf/metrics/qso_number_counts_metric.py +++ b/rubin_sim/maf/metrics/qso_number_counts_metric.py @@ -12,13 +12,16 @@ class QSONumberCountsMetric(BaseMetric): - """ - Calculate the number of quasars expected with SNR>=5 according to the Shen et al. (2020) QLF - - model A in the redshift range zmin < z < zmax. The 5 sigma depths are obtained using the ExgalM5 metric. + """Calculate the number of quasars expected with SNR>=5 + according to the Shen et al. (2020) QLF - model A in the redshift + range zmin < z < zmax. + + The 5 sigma depths are obtained using the ExgalM5 metric. Only quasars fainter than the saturation magnitude are counted. - By default, zmin is 0.3 and zmax is the minimum between 6.7 and the redshift at which the Lyman break - matches the effective wavelength of the band. For bands izy, zmax is 6.7. This default choice is to + By default, zmin is 0.3 and zmax is the minimum between 6.7 and the + redshift at which the Lyman break matches the effective wavelength + of the band. For bands izy, zmax is 6.7. This default choice is to match Table 10.2 for i-band quasar counts in the LSST Science book. """ @@ -47,7 +50,8 @@ def __init__( "y": 971.0, } - # Dust Extinction limit. Regions with larger extinction and dropped from the counting. + # Dust Extinction limit. + # Regions with larger extinction and dropped from the counting. self.extinction_cut = extinction_cut # Save the filter information. @@ -55,9 +59,11 @@ def __init__( self.lsst_filter = lsst_filter # Save zmin and zmax, or set zmax to the default value. - # The default zmax is the lower number between 6.7 and the redshift at which the - # Lyman break (91.2nm) hits the effective wavelength of the filter. - # Note that this means that for i, z and y the default value for zmax is 6.7 + # The default zmax is the lower number between 6.7 and the + # redshift at which the Lyman break (91.2nm) hits the + # effective wavelength of the filter. + # Note that this means that for i, z and y, + # the default value for zmax is 6.7 self.zmin = zmin if zmax is None: zmax = np.min([6.7, self.effwavelen[self.lsst_filter] / 91.2 - 1.0]) @@ -71,8 +77,9 @@ def __init__( self.qlf_model = qlf_model self.sed_model = sed_model - # Read the long tables, which the number of quasars expected for a given band, - # qlf_module and qlf_model in a range of redshifts and magnitudes. + # Read the long tables, which the number of quasars expected + # for a given band, qlf_module and qlf_model in a range of + # redshifts and magnitudes. table_name = "Long_Table.LSST{0}.{1}.{2}.{3}.txt".format( self.lsst_filter, self.qlf_module, self.qlf_model, self.sed_model ) @@ -90,7 +97,8 @@ def __init__( c_mz_data = np.cumsum(c_mz_data, axis=1) # Create a 2D interpolation object for the long table. - # self.nqso_cumulative = interpolate.interp2d(zs[:-1], mags[:-1], #c_mz_data[:-1, :-1], kind="cubic") + # self.nqso_cumulative = interpolate.interp2d(zs[:-1], mags[:-1], + # #c_mz_data[:-1, :-1], kind="cubic") self.nqso_cumulative_aux = interpolate.RectBivariateSpline( zs[:-1], mags[:-1], c_mz_data[:-1, :-1].T, kx=3, ky=3 ) diff --git a/rubin_sim/maf/metrics/scaling_metrics.py b/rubin_sim/maf/metrics/scaling_metrics.py index cb25a1ca5..1589000f7 100644 --- a/rubin_sim/maf/metrics/scaling_metrics.py +++ b/rubin_sim/maf/metrics/scaling_metrics.py @@ -14,13 +14,13 @@ class NgalScaleMetric(BaseMetric): Parameters ---------- - a_max : float (0.2) + a_max : `float` The maximum dust extinction to allow. Anything with higher dust extinction is considered to have zero usable galaxies. - m5min : float (26) + m5min : `float` The minimum coadded 5-sigma depth to allow. Anything less is considered to have zero usable galaxies. - filter : str ("i") + filter : `str` The filter to use. Any visits in other filters are ignored. """ @@ -55,7 +55,8 @@ def __init__( self.ax1 = dust_properties.ax1 def run(self, data_slice, slice_point): - # I'm a little confused why there's a dust cut and an M5 cut, but whatever + # I'm a little confused why there's a dust cut and an M5 cut, + # but whatever a_x = self.ax1[data_slice[self.filter_col][0]] * slice_point["ebv"] if a_x > self.a_max: return 0 @@ -80,16 +81,18 @@ class NlcPointsMetric(BaseMetric): Parameters ---------- - ndpmin : int (10) + ndpmin : `int` The number of points to demand on a lightcurve in a single filter to have that light curve qualify. - mags : float (21) + mags : `float` The magnitude of our fiducial object (maybe make it a dict in the future to support arbitrary colors). - maps : list of map objects (None) - List of stellar density maps to use. Default of None loads Trilegal maps. - nside : int (128) - The nside is needed to make sure the loaded maps match the slicer nside. + maps : `list` [`~rubin_sim.maf.map`] or None + List of stellar density maps to use. + Default of None loads Trilegal maps. + nside : `int` + The nside is needed to make sure the loaded maps + match the slicer nside. """ def __init__( diff --git a/rubin_sim/maf/metrics/season_metrics.py b/rubin_sim/maf/metrics/season_metrics.py index f386990c4..229d0ad56 100644 --- a/rubin_sim/maf/metrics/season_metrics.py +++ b/rubin_sim/maf/metrics/season_metrics.py @@ -1,5 +1,6 @@ -"""A group of metrics that work together to evaluate season characteristics (length, number, etc). -In addition, these supports the time delay metric calculation for strong lensing. +"""A group of metrics that work together to evaluate season characteristics +(length, number, etc). +In addition, these support the time delay metric for strong lensing. """ __all__ = ( "find_season_edges", @@ -22,13 +23,13 @@ def find_season_edges(seasons): Parameters ---------- - seasons: np.ndarray + seasons : `np.ndarray`, (N,) Seasons, such as calculated by calc_season. Note that seasons should be sorted!! Returns ------- - np.ndarray, np.ndarray + first, last : `np.ndarray`, (N,), `np.ndarray`, (N,) The indexes of the first and last date in the season. """ int_seasons = np.floor(seasons) @@ -41,17 +42,24 @@ def find_season_edges(seasons): class SeasonLengthMetric(BaseMetric): - """ - Calculate the length of LSST seasons, in days. + """Calculate the length of LSST seasons, in days. Parameters ---------- - min_exp_time: float, optional - Minimum visit exposure time to count for a 'visit', in seconds. Default 20. + min_exp_time : `float`, optional + Minimum visit exposure time to count for a 'visit', in seconds. + Default 20. reduce_func : function, optional - Function that can operate on array-like structures. Typically numpy function. - This reduces the season length in each season from 10 separate values to a single value. + Function that can operate on array-like structures. + Typically numpy function. + This reduces the season length in each season from 10 separate + values to a single value. Default np.median. + + Returns + ------- + seasonlength : `float` + The (reduceFunc) of the length of each season, in days. """ def __init__( @@ -73,28 +81,13 @@ def __init__( ) def run(self, data_slice, slice_point): - """Calculate the (reduceFunc) of the length of each season. - Uses the slice_point RA/Dec to calculate the position in question, then uses the times of the visits - to assign them into seasons (based on where the sun is relative to the slice_point RA). - - Parameters - ---------- - data_slice : numpy.array - Numpy structured array containing the data related to the visits provided by the slicer. - slice_point : dict - Dictionary containing information about the slice_point currently active in the slicer. - - Returns - ------- - float - The (reduceFunc) of the length of each season, in days. - """ # Order data Slice/times and exclude visits which are too short. long = np.where(data_slice[self.exp_time_col] > self.min_exp_time) if len(long[0]) == 0: return self.badval data = np.sort(data_slice[long], order=self.mjd_col) - # SlicePoints ra/dec are always in radians - convert to degrees to calculate season + # SlicePoints ra/dec are always in radians - + # convert to degrees to calculate season seasons = calc_season(np.degrees(slice_point["ra"]), data[self.mjd_col]) first_of_season, last_of_season = find_season_edges(seasons) seasonlengths = data[self.mjd_col][last_of_season] - data[self.mjd_col][first_of_season] @@ -103,7 +96,8 @@ def run(self, data_slice, slice_point): class CampaignLengthMetric(BaseMetric): - """Calculate the number of seasons (roughly, years) a pointing is observed for. + """Calculate the number of seasons (roughly, years) a pointing is observed. + This corresponds to the 'campaign length' for lensed quasar time delays. """ @@ -129,7 +123,9 @@ def run(self, data_slice, slice_point): class MeanCampaignFrequencyMetric(BaseMetric): - """Calculate the mean separation between nights, within a season - then the mean over the campaign. + """Calculate the mean separation between nights, within a season - + then the mean over the campaign. + Calculate per season, to avoid any influence from season gaps. """ @@ -154,7 +150,8 @@ def run(self, data_slice, slice_point): if len(long[0]) == 0: return self.badval data = np.sort(data_slice[long], order=self.mjd_col) - # SlicePoints ra/dec are always in radians - convert to degrees to calculate season + # SlicePoints ra/dec are always in radians - + # convert to degrees to calculate season seasons = calc_season(np.degrees(slice_point["ra"]), data[self.mjd_col]) first_of_season, last_of_season = find_season_edges(seasons) season_means = np.zeros(len(first_of_season), float) @@ -168,49 +165,54 @@ def run(self, data_slice, slice_point): class TdcMetric(BaseMetric): - """Calculate the Time Delay Challenge metric, as described in Liao et al 2015 - (https://arxiv.org/pdf/1409.1254.pdf). + """Calculate the Time Delay Challenge metric, + as described in Liao et al 2015 (https://arxiv.org/pdf/1409.1254.pdf). - This combines the MeanCampaignFrequency/MeanNightSeparation, the SeasonLength, and the CampaignLength + This combines the MeanCampaignFrequency/MeanNightSeparation, + the SeasonLength, and the CampaignLength metrics above, but rewritten to calculate season information only once. cad_norm = in units of days sea_norm = in units of months camp_norm = in units of years - This metric also adds a requirement to achieve limiting magnitudes after galactic dust extinction, - in various bandpasses, in order to exclude visits which are not useful for detecting quasars - (due to being short or having high sky brightness, etc.) and to reject regions with - high galactic dust extinction. + This metric also adds a requirement to achieve limiting magnitudes + after galactic dust extinction, in various bandpasses, + in order to exclude visits which are not useful for detecting quasars + (due to being short or having high sky brightness, etc.) and to + reject regions with high galactic dust extinction. Parameters ---------- - mjd_col: str, optional + mjd_col : `str`, optional Column name for mjd. Default observationStartMJD. - night_col: str, optional + night_col : `str`, optional Column name for night. Default night. - filter_col: str, optional + filter_col : `str`, optional Column name for filter. Default filter. - m5_col: str, optional + m5_col : `str`, optional Column name for five-sigma depth. Default fiveSigmaDepth. - mag_cuts: dict, optional - Dictionary with filtername:mag limit (after dust extinction). Default None in kwarg. - Defaults set within metric: {'u': 22.7, 'g': 24.1, 'r': 23.7, 'i': 23.1, 'z': 22.2, 'y': 21.4} - metricName: str, optional + mag_cuts : `dict`, optional + Dictionary with filtername:mag limit (after dust extinction). + Default None in kwarg. + Defaults set within metric: + {'u': 22.7, 'g': 24.1, 'r': 23.7, 'i': 23.1, 'z': 22.2, 'y': 21.4} + metricName : `str`, optional Metric Name. Default TDC. - cad_norm: float, optional + cad_norm : `float`, optional Cadence normalization constant, in units of days. Default 3. - sea_norm: float, optional + sea_norm : `float`, optional Season normalization constant, in units of months. Default 4. - camp_norm: float, optional + camp_norm : `float`, optional Campaign length normalization constant, in units of years. Default 5. - badval: float, optional + badval : `float`, optional Return this value instead of the dictionary for bad points. Returns ------- - dictionary - Dictionary of values for {'rate', 'precision', 'accuracy'} at this point in the sky. + TDCmetrics : `dict` + Dictionary of values for {'rate', 'precision', 'accuracy'} + at this point in the sky. """ def __init__( @@ -250,7 +252,8 @@ def __init__( raise Exception("mag_cuts should be a dictionary") # Set up dust map requirement maps = ["DustMap"] - # Set the default wavelength limits for the lsst filters. These are approximately correct. + # Set the default wavelength limits for the lsst filters. + # These are approximately correct. dust_properties = DustValues() self.ax1 = dust_properties.ax1 super().__init__( @@ -275,7 +278,8 @@ def run(self, data_slice, slice_point): if len(idxs[0]) == 0: return self.badval data = np.sort(data_slice[idxs], order=self.mjd_col) - # SlicePoints ra/dec are always in radians - convert to degrees to calculate season + # SlicePoints ra/dec are always in radians - + # convert to degrees to calculate season seasons = calc_season(np.degrees(slice_point["ra"]), data[self.mjd_col]) int_seasons = np.floor(seasons) first_of_season, last_of_season = find_season_edges(seasons) diff --git a/rubin_sim/maf/metrics/simple_metrics.py b/rubin_sim/maf/metrics/simple_metrics.py index 134f12cad..2621f51e4 100644 --- a/rubin_sim/maf/metrics/simple_metrics.py +++ b/rubin_sim/maf/metrics/simple_metrics.py @@ -290,7 +290,9 @@ def run(self, data_slice, slice_point=None): class MaxPercentMetric(BaseMetric): - """Return the percent of data which matches the maximum value of the data.""" + """Return the percent of data which matches the maximum value + of the data. + """ def run(self, data_slice, slice_point=None): n_max = np.size(np.where(data_slice[self.colname] == np.max(data_slice[self.colname]))[0]) @@ -324,8 +326,9 @@ class FracAboveMetric(BaseMetric): """Find the fraction of data values above a given `cutoff`.""" def __init__(self, col=None, cutoff=0.5, scale=1, metric_name=None, **kwargs): - # Col could just get passed in bundle with kwargs, but by explicitly pulling it out - # first, we support use cases where class instantiated without explicit 'col='). + # Col could just get passed in bundle with kwargs, + # by explicitly pulling it out first, we support use cases where + # class instantiated without explicit 'col='). if metric_name is None: metric_name = "FracAbove %.2f in %s" % (cutoff, col) super(FracAboveMetric, self).__init__(col, metric_name=metric_name, **kwargs) diff --git a/rubin_sim/maf/metrics/slew_metrics.py b/rubin_sim/maf/metrics/slew_metrics.py deleted file mode 100644 index 7484079f3..000000000 --- a/rubin_sim/maf/metrics/slew_metrics.py +++ /dev/null @@ -1,71 +0,0 @@ -__all__ = ("SlewContributionMetric", "AveSlewFracMetric") - -import numpy as np - -from .base_metric import BaseMetric - -# Metrics for dealing with things from the SlewActivities table - - -class SlewContributionMetric(BaseMetric): - def __init__( - self, col="actDelay", activity=None, active_col="activity", in_crit_col="inCriticalPath", **kwargs - ): - """ - Return the average time, multiplied by fraction of slew -- - considering critical path activities only. - """ - self.col = col - self.in_crit_col = in_crit_col - col = [col, in_crit_col] - col.append(active_col) - self.active_col = active_col - self.activity = activity - super(SlewContributionMetric, self).__init__(col=col, **kwargs) - self.comment = "Average time for %s activity (in seconds) when in the critical path, " % (activity) - self.comment += "multiplied by the percent of total slews in the critical path." - - def run(self, data_slice, slice_point=None): - # Activities of this type, in critical path. - good_in_crit = np.where( - (data_slice[self.active_col] == self.activity) & (data_slice[self.in_crit_col] == "True") - )[0] - if len(good_in_crit) == 0: - result = 0.0 - else: - # All activities in critical path. - in_crit = np.where((data_slice[self.in_crit_col] == "True"))[0] - # Calculate fraction of total in-critical-path slew activities that this activity represents. - result = np.sum(data_slice[self.col][good_in_crit]) / np.sum(data_slice[self.col][in_crit]) - # and multiply by the mean time required by this activity. - result *= np.mean(data_slice[self.col][good_in_crit]) - return result - - -class AveSlewFracMetric(BaseMetric): - def __init__( - self, col="actDelay", activity=None, active_col="activity", id_col="SlewHistory_slewCount", **kwargs - ): - """ - Return the average time multiplied by fraction of slews. - """ - self.col = col - self.id_col = id_col - col = [col, id_col] - col.append(active_col) - self.active_col = active_col - self.activity = activity - super(AveSlewFracMetric, self).__init__(col=col, **kwargs) - self.comment = "Average time for %s activity (in seconds), multiplied by percent of total slews." % ( - activity - ) - - def run(self, data_slice, slice_point=None): - good = np.where(data_slice[self.active_col] == self.activity)[0] - if len(good) == 0: - result = 0.0 - else: - result = np.mean(data_slice[self.col][good]) - nslews = np.size(np.unique(data_slice[self.id_col])) - result = result * np.size(good) / np.float(nslews) - return result diff --git a/rubin_sim/maf/metrics/sn_cadence_metric.py b/rubin_sim/maf/metrics/sn_cadence_metric.py index 580d91ca1..201ae2cd7 100644 --- a/rubin_sim/maf/metrics/sn_cadence_metric.py +++ b/rubin_sim/maf/metrics/sn_cadence_metric.py @@ -7,16 +7,17 @@ class SNCadenceMetric(metrics.BaseMetric): """ - Metric to estimate the redshift limit for faint supernovae (x1,color) = (-2.0,0.2) + Metric to estimate the redshift limit for faint supernovae + (x1,color) = (-2.0,0.2) Parameters ---------- - list : str, optional + list : `str`, optional Name of the columns used to estimate the metric - coadd : bool, optional + coadd : `bool`, optional to make "coaddition" per night (uses snStacker) Default True - lim_sn : class, optional + lim_sn : `class`, optional Reference data used to estimate redshift values (interpolation) """ diff --git a/rubin_sim/maf/metrics/sn_n_sn_metric.py b/rubin_sim/maf/metrics/sn_n_sn_metric.py index c7103e456..c7edf7244 100644 --- a/rubin_sim/maf/metrics/sn_n_sn_metric.py +++ b/rubin_sim/maf/metrics/sn_n_sn_metric.py @@ -69,7 +69,8 @@ class SNNSNMetric(BaseMetric): dust : `bool`, opt Apply dust extinction to visit depth values (default False) hard_dust_cut : `float`, opt - If set, cut any point on the sky that has an ebv extinction higher than the hard_dust_cut value. + If set, cut any point on the sky that has an ebv extinction + higher than the hard_dust_cut value. Default 0.25 """ @@ -221,14 +222,14 @@ def run(self, data_slice, slice_point): Parameters -------------- - data_slice: `np.array` + data_slice : `np.ndarray` Observations to process (scheduler simulations) - slice_point: `bool`, opt + slice_point : `bool`, opt Information about the location on the sky from the slicer Returns ------- - metricVal : `np.recarray` + metricVal : `np.ndarray` ['n_sn', 'zlim'] at this point on the sky """ # Hard dust cut @@ -264,7 +265,8 @@ def run(self, data_slice, slice_point): if len(data_slice) <= self.n_aft + self.n_bef: return self.badval - # get season information (seasons calculated by gaps, not by place on sky) + # get season information (seasons calculated by gaps, + # not by place on sky) data_slice = self.getseason(data_slice, mjd_col=self.mjd_col) # get redshift values per season @@ -326,17 +328,17 @@ def season_length(self, seasons, data_slice, zseason): Method to estimate season lengths vs z Parameters - --------------- - seasons : list(int) + ----------- + seasons : `list` [`int`] list of seasons to process - data_slice: numpy array + data_slice : `np.ndarray`, (N,)` array of observations Returns - ----------- - seasons : list(int) + -------- + seasons : `list` [`int`] list of seasons to process - dur_z : pandas df + dur_z : `pandas.DataFrame` season lengths vs z """ # if seasons = -1: process the seasons seen in data @@ -368,15 +370,15 @@ def get_season_info(self, dfa, zseason, min_duration=60.0): Parameters -------------- - dfa: pandas df + dfa : pandas df dat to process - zseason: pandas df + zseason : pandas df redshift infos per season - min_duration: float, opt + min_duration : `float`, opt min season length to be accepted (default: 60 days) Returns - ---------- + -------- pandas df with season length infos """ @@ -406,13 +408,13 @@ def step_lc(self, obs, gen_par, x1=-2.0, color=0.2): Parameters --------------- - obs: array + obs : array observations - gen_par: array + gen_par : array simulation parameters - x1: float, opt + x1 : `float`, opt stretch value (default: -2.0) - color: float, opt + color : `float`, opt color value (default: 0.2) Returns @@ -437,9 +439,6 @@ def step_efficiencies(self, lc): pandas df with efficiencies """ - # sn_effis = lc.groupby(['healpixID', 'season', 'z', 'x1', 'color', 'sntype']).apply( - # lambda x: self.sn_effi(x)).reset_index() - sn_effis = ( lc.groupby(["season", "z", "x1", "color", "sntype"]) .apply(lambda x: self.sn_effi(x)) @@ -476,14 +475,14 @@ def step_nsn(self, sn_effis, dur_z): Method to estimate the number of supernovae from efficiencies Parameters - --------------- - sn_effis: pandas df + ---------- + sn_effis : pandas df data with efficiencies of observation - dur_z: array + dur_z : array array of season length Returns - ---------- + ------- initial sn_effis appended with a set of infos (duration, nsn) """ @@ -501,8 +500,8 @@ def season_info(self, grp, min_duration): Parameters -------------- - grp: pandas df group - min_duration: float + grp : pandas df group + min_duration : `float` minimal duration for a season to be considered Returns @@ -543,9 +542,9 @@ def duration_z(self, grp, min_duration=60.0): Parameters -------------- - grp: pandas df group + grp : pandas df group data to process: season infos - min_duration: float, opt + min_duration : `float`, opt min season length for a season to be processed (deafult: 60 days) Returns @@ -581,11 +580,11 @@ def calc_daymax(self, grp, daymax_step): Parameters -------------- grp: group (pandas df sense) - group of data to process with the following cols: - t0_min: T0 min value (per season) - t0_max: T0 max value (per season) - daymax_step: float - step for T0 simulation + group of data to process with the following cols: + t0_min: T0 min value (per season) + t0_max: T0 max value (per season) + daymax_step: `float` + step for T0 simulation Returns ---------- @@ -615,13 +614,13 @@ def gen_lc(self, grp, gen_par_orig, x1, color): Parameters --------------- - grp: pandas group + grp : pandas group observations to process - gen_par_orig: pandas df + gen_par_orig : pandas df simulation parameters - x1: float + x1 : `float` SN stretch - color: float + color : `float` SN color Returns @@ -656,7 +655,7 @@ def sn_effi(self, lc): Parameters --------------- - lc: pandas grp + lc : pandas grp light curve Returns @@ -721,13 +720,13 @@ def get_sum(self, lcarr, varname, nvals, flag): Parameters -------------- - lcarr: numpy array + lcarr : numpy array data to process - varname: str + varname : `str` col to process in lcarr - nvals: int + nvals : `int` dimension for tiling - flag: array(bool) + flag : array(bool) flag to apply Returns @@ -748,11 +747,11 @@ def get_epochs(self, nights, flag, flagph): Parameters --------------- - nights: array + nights : array night number array - flag: array(bool) + flag : array(bool) flag to apply - flagph: array(bool) + flagph : array(bool) flag to apply Returns @@ -1316,7 +1315,8 @@ def get_nsn(self, effi, durinterp_z, zmin, zmax, zstep): def check_dur_z(self, dur_z, nmin=2): """ " - Method to remove seasons with a poor redshift range due to too low season length + Method to remove seasons with a poor redshift range due + to too low season length Parameters ---------------- diff --git a/rubin_sim/maf/plots/spatial_plotters.py b/rubin_sim/maf/plots/spatial_plotters.py index b0f15ef8d..02fae6b07 100644 --- a/rubin_sim/maf/plots/spatial_plotters.py +++ b/rubin_sim/maf/plots/spatial_plotters.py @@ -145,16 +145,18 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None): """ Parameters ---------- - metric_value : numpy.ma.MaskedArray - slicer : rubin_sim.maf.slicers.HealpixSlicer - user_plot_dict: dict - Dictionary of plot parameters set by user (overrides default values). - fignum : int - Matplotlib figure number to use (default = None, starts new figure). + metric_value : `numpy.ma.MaskedArray` + slicer : `rubin_sim.maf.slicers.HealpixSlicer` + user_plot_dict: `dict` + Dictionary of plot parameters set by user + (overrides default values). + fignum : `int` + Matplotlib figure number to use + (default = None, starts new figure). Returns ------- - int + fignum : `int` Matplotlib figure number used to create the plot. """ # Override the default plotting parameters with user specified values. @@ -248,7 +250,8 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None): visufunc_params.update(self.healpy_visufunc_params) self.healpy_visufunc(metric_value.filled(badval), **visufunc_params) - # Add colorbar (not using healpy default colorbar because we want more tickmarks). + # Add colorbar + # (not using healpy default colorbar because we want more tickmarks). self.ax = plt.gca() im = self.ax.get_images()[0] @@ -259,12 +262,12 @@ def __call__(self, metric_value_in, slicer, user_plot_dict, fignum=None): # Add label. if plot_dict["label"] is not None: plt.figtext(0.8, 0.8, "%s" % (plot_dict["label"])) - # Make a color bar. Supress silly colorbar warnings. + # Make a color bar. Suppress excessive colorbar warnings. with warnings.catch_warnings(): warnings.simplefilter("ignore") # The vertical colorbar is primarily aimed at the movie # but may be useful for other purposes - if plot_dict["extend"] is not "neither": + if plot_dict["extend"] != "neither": extendrect = False else: extendrect = True @@ -315,8 +318,8 @@ def __init__(self): self.default_plot_dict.update({"maxl": None, "removeDipole": True, "linestyle": "-"}) def __call__(self, metric_value, slicer, user_plot_dict, fignum=None): - """ - Generate and plot the power spectrum of metric_value (calculated on a healpix grid). + """Generate and plot the power spectrum of metric_values + (for metrics calculated on a healpix grid). """ if "Healpix" not in slicer.slicer_name: raise ValueError("HealpixPowerSpectrum for use with healpix metricBundles.") @@ -361,7 +364,8 @@ def __call__(self, metric_value, slicer, user_plot_dict, fignum=None): plt.tick_params(axis="y", labelsize=plot_dict["labelsize"]) if plot_dict["title"] is not None: plt.title(plot_dict["title"]) - # Return figure number (so we can reuse/add onto/save this figure if desired). + # Return figure number + # (so we can reuse/add onto/save this figure if desired). return fig.number @@ -384,8 +388,7 @@ def __init__(self): self.base_hist = BaseHistogram() def __call__(self, metric_value, slicer, user_plot_dict, fignum=None): - """ - Histogram metric_value for all healpix points. + """Histogram metric_value for all healpix points. """ if "Healpix" not in slicer.slicer_name: raise ValueError("HealpixHistogram is for use with healpix slicer.") diff --git a/rubin_sim/maf/stackers/coord_stackers.py b/rubin_sim/maf/stackers/coord_stackers.py index 1605184ff..b3cf93ac6 100644 --- a/rubin_sim/maf/stackers/coord_stackers.py +++ b/rubin_sim/maf/stackers/coord_stackers.py @@ -13,28 +13,28 @@ def ra_dec2_alt_az(ra, dec, lat, lon, mjd, altonly=False): """Convert RA/Dec (and telescope site lat/lon) to alt/az. - This uses simple equations and ignores aberation, precession, nutation, etc. + This uses simple equations and ignores aberation, precession, nutation. Parameters ---------- - ra : array_like + ra : `np.ndarray`, (N,) RA, in radians. - dec : array_like + dec : `np.ndarray`, (N,) Dec, in radians. Must be same length as `ra`. - lat : float + lat : `float` Latitude of the observatory in radians. - lon : float + lon : `float` Longitude of the observatory in radians. - mjd : float + mjd : `float` Modified Julian Date. - altonly : bool, optional + altonly : `bool`, optional Calculate altitude only. Returns ------- - alt : numpy.array + alt : `np.ndarray`, (N,) Altitude, same length as `ra` and `dec`. Radians. - az : numpy.array + az : `np.ndarray`, (N,) Azimuth, same length as `ra` and `dec`. Radians. """ lmst = calc_lmst(mjd, lon) diff --git a/rubin_sim/maf/utils/generate_fov_map.py b/rubin_sim/maf/utils/generate_fov_map.py index eabea3b11..36afda59b 100644 --- a/rubin_sim/maf/utils/generate_fov_map.py +++ b/rubin_sim/maf/utils/generate_fov_map.py @@ -2,6 +2,7 @@ from rubin_scheduler.utils import gnomonic_project_tosky, gnomonic_project_toxy # Use the main stack to make a rough array. +# This code needs an update to work without lsst.sims. if __name__ == "__main__": import lsst.sims.utils as simsUtils diff --git a/rubin_sim/maf/web/maf_run_results.py b/rubin_sim/maf/web/maf_run_results.py index d4fd42f5b..39fd6fcab 100644 --- a/rubin_sim/maf/web/maf_run_results.py +++ b/rubin_sim/maf/web/maf_run_results.py @@ -177,7 +177,6 @@ def metric_ids_to_metrics(self, metric_ids, metrics=None): """ if metrics is None: metrics = self.metrics - # this should be faster with pandas (and self.metrics.query('metric_id in @metric_ids')) metrics = metrics[np.in1d(metrics["metric_id"], metric_ids)] return metrics @@ -449,7 +448,8 @@ def plot_dict(self, plots=None): Given an array of plots (for a single metric usually). Returns an ordered dict with 'plot_type' for interfacing with jinja2 templates. - plot_dict == {'SkyMap': {'plot_file': [], 'thumb_file', []}, 'Histogram': {}..} + plot_dict == + {'SkyMap': {'plot_file': [], 'thumb_file', []}, 'Histogram': {}..} If no plot of a particular type, the plot_file and thumb_file are empty lists.