#!/usr/bin/env python3
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi
# Date: 2023-07-12 09:32:38
"""General utilities."""
import logging
import os
import astropy.units as u
import fitsio
import numpy as np
from astropy import constants
from astropy.coordinates import SkyCoord
from astropy.cosmology import FlatLambdaCDM
from astropy.time import Time
from scipy.stats import binned_statistic
from lsst_inaf_agile import zou2024
from lsst_inaf_agile.igm_absorption import my_get_IGM_absorption
logging.basicConfig(format="%(asctime)s - %(name)s - %(levelname)s - %(message)s", level=logging.INFO)
[docs]
logger = logging.getLogger(__name__)
[docs]
def flux_to_mag(flux):
"""
Convert uJy flux to AB magnitude.
Parameters
----------
flux: float
Input flux in microjanskies.
Examples
--------
>>> import numpy as np
>>> "%.2f" % flux_to_mag(1.0)
'23.90'
>>> "%.2f" % flux_to_mag(10.)
'21.40'
>>> flux_to_mag(0.)
array(nan)
"""
return np.where(np.atleast_1d(flux) > 0, -2.5 * np.ma.log10(flux * 1e-6 / 3631), np.nan).squeeze()
[docs]
def mag_to_flux(mag):
"""
Convert AB magnitude to uJy flux.
Parameters
----------
mag: float
Input AB magnitude.
Examples
--------
>>> import numpy as np
>>> "%.2f" % mag_to_flux(23.90)
'1.00'
>>> "%.2f" % mag_to_flux(21.40)
'10.00'
>>> mag_to_flux(np.inf)
0.0
>>> "%.2f" % (mag_to_flux(0) / 1e6)
'3631.00'
"""
return 3631 * 1e6 * 10 ** (mag / -2.5)
[docs]
def mag_sum(mag):
"""
Sum together magnitudes 'mag', return the combined magnitude.
Parameters
----------
mag: float
Magnitude(s) to sum together.
Examples
--------
>>> import numpy as np
>>> mag_sum(0.0)
array(-0.)
>>> mag_sum([0.0, 2.5])
array(-0.10348171)
>>> mag_sum([0.0, 1.0, 2.0])
array(-0.48044012)
"""
flux = np.sum(mag_to_flux(np.array(mag)))
return flux_to_mag(flux)
[docs]
def get_volume(
zmin: float,
zmax: float,
area_deg2: float = 41252.96124941928,
H0: float = 70.0,
Om0: float = 0.30,
Tcmb0: float = 2.73,
):
"""
Return the flat LambdaCDM comoving volume in Mpc for a redshift shell.
Parameters
----------
zmin: float
Minimum redshift.
zmax: float
Maximum redshift.
area_deg2: float
Sky area in square degrees.
H0: float
Present day Hubble parameter in km/s/Mpc. Use H0=100 (i.e. h=1) to
return the comoving volume in (Mpc/h)^3
Om0: float
Present day dimensionless matter density parameter.
Tcmb0: float
Present day CMB temperature in Kelvin.
Returns
-------
comoving_volume: float
The comoving volume segment in Mpc3 corresponding to the input
arguments.
Examples
--------
>>> # Following values only illustrate the usage.
>>> # Accuracy of the volume estimation is set by AstroPy's implementation.
>>> get_volume(0.0, 1.0, 1.0)
np.float64(3660715.356254536)
>>> get_volume(1.0, 2.0, 1.0)
np.float64(10444274.253266422)
>>> get_volume(0.0, 1.0, 10.0)
np.float64(36607153.56254536)
>>> get_volume(0.0, 1.0, 1.0, Om0=0.45)
np.float64(2887067.685807227)
>>> get_volume(0.0, 1.0, 1.0, Tcmb0=0.00)
np.float64(3661728.0492256973)
"""
cosmo = FlatLambdaCDM(H0=H0, Om0=Om0, Tcmb0=Tcmb0)
volume = (cosmo.comoving_volume(zmax) - cosmo.comoving_volume(zmin)).value
ret = volume * area_deg2 / (4 * np.pi * u.sr.to(u.deg**2))
return ret
[docs]
def get_chisq_nu(y, y_model, sigma):
"""
Return the chisq value of the measurement y.
Input arguments are converted to numpy arrays before calculation.
Parameters
----------
y: float
Measured values.
y_model: float
Model values.
sigma: float
Model errors.
Returns
-------
chisq_nu: float
Estimated reduced chisq value. Degrees of freedom are assumed to be
len(y)-1.
Raises
------
ValueError
If number of good (finite) data points is less than two.
Examples
--------
>>> import numpy as np
>>> get_chisq_nu(1, 1, 1)
Traceback (most recent call last):
...
ValueError: Number of good data point is less than two.
>>> get_chisq_nu([1, 2, 3], [3, 2, 1], [1, 3, 2])
np.float64(2.5)
>>> get_chisq_nu([1, 2, np.nan], [3, 2, np.nan], [1, 3, np.nan])
np.float64(4.0)
"""
y = np.array(y)
y_model = np.array(y_model)
sigma = np.array(sigma)
diff = np.ma.true_divide(y - y_model, sigma)
is_good = np.isfinite(y) & np.isfinite(y_model) & np.isfinite(sigma)
if np.sum(is_good) < 2:
raise ValueError("Number of good data point is less than two.")
return np.ma.true_divide(np.ma.sum(diff[is_good] ** 2), is_good.sum() - 1)
[docs]
def get_key_function(bins, x, values=None, nmin=30, *args, **kwargs):
"""
Return the "key" function i.e. number of objects at interval in 'key'.
Key is e.g the stellar mass or the X-ray luminosity. The function is
returned in units of 1/Mpc3/dex for the given cosmology.
Parameters
----------
bins: list[float]
Edges of the bins.
x: list[float]
Values to be binned.
values: list[float] or None
Weights to 'x'. Can be None in which case no weighting is done.
nmin: int
Minimum number of counts per bin to be considered 'valid'. Bins with
less than 'nmin' counts are masked to negative values.
args, kwargs:
Additional arguments forwarded to the function 'get_volume'.
Returns
-------
x, dx, y, dy:
Bin centers (x), Bin width (dx), Counts (y), Delta counts (dy). Error
on the counts is assumed to be Poissonian i.e. sqrt(Ncounts).
Examples
--------
>>> import numpy as np
>>> bins = np.array([0, 1])
>>> # Less than nmin counts returns negative values
>>> x, dx, y, dy = get_key_function(bins, np.array([0]), zmin=0.00, zmax=0.10)
>>> assert np.all(y < 0)
>>> # More than nmin counts returns positive values
>>> x, dx, y, dy = get_key_function(bins, np.array([0] * 31), zmin=0.00, zmax=0.10)
>>> y[0]
np.float64(1.0100431160959446e-07)
>>> # Doubling the values doubles the returned function
>>> x, dx, y, dy = get_key_function(bins, np.array([0] * 31 * 2), zmin=0.00, zmax=0.10)
>>> y[0]
np.float64(2.0200862321918891e-07)
>>> # Weighting modifies the output data
>>> x, dx, y, dy = get_key_function(
... bins, np.array([0] * 31 * 2), values=np.array([1] * 31 + [0] * 31), zmin=0.00, zmax=0.10
... )
>>> y[0]
np.float64(1.0100431160959446e-07)
>>> # Empty array returns a zero
>>> x, dx, y, dy = get_key_function(bins, np.array([]), zmin=0.00, zmax=0.10)
>>> y[0]
np.float64(0.0)
"""
# The binning
dbins = np.diff(bins)
assert np.allclose(dbins[0], dbins[1:])
centers = bins[:-1] + dbins / 2
# Short-circuit for an empty array
if not x.size:
return centers, dbins, np.zeros_like(centers), np.zeros_like(centers)
if values is None:
values = np.ones_like(x)
# Calculate the weight by using the area and the volume
counts = binned_statistic(x=x, values=values, statistic="sum", bins=bins)[0]
counts = np.where(counts > nmin, counts, -99)
div = get_volume(*args, **kwargs) * dbins[0]
function = np.ma.true_divide(counts, div)
# NOTE: assume poissonian errors on the counts, no cosmological variance etc.
dfunction = np.ma.sqrt(counts) / div
return centers, dbins / 2, function, dfunction
[docs]
def egg_band_to_index(egg: dict, band: str) -> int:
"""
Convert EGG band name to an index.
Parameters
----------
egg: dict
Dictionary-like EGG-dataset. Can be simple output from reading an EGG
FITS file.
band: str
Name of the band e.g. 'lsst-r'.
Returns
-------
idx: int
Index of the band in the EGG catalog.
Raises
------
ValueError
If EGG does not contain the band.
Examples
--------
>>> # Generate a mock EGG catalog -- note the shape [1, 3]
>>> egg = {"BANDS": [["lsst-r", "lsst-g", "lsst-i"]]}
>>> egg_band_to_index(egg, "lsst-r")
0
>>> egg_band_to_index(egg, "lsst-g")
1
>>> egg_band_to_index(egg, "lsst-i")
2
>>> egg_band_to_index(egg, "non-existing-band")
Traceback (most recent call last):
...
ValueError: 'non-existing-band' is not in list
"""
bands = [b.strip() for b in egg["BANDS"][0]]
return bands.index(band)
[docs]
def get_ra_dec(ra0, dec0, pm_ra_cosdec, pm_dec, mjd, mjd0=51544.5):
"""
Get ra, dec for current epoch modified by the proper motion.
Parameters
----------
ra0: float
Right ascension at mjd0.
dec0: float
Declination at mjd0.
pm_ra_cosdec: float
Proper motion in (right ascension) * cos(declination) in mas/yr.
pm_dec: float
Proper motion in declination in mas/yr.
mjd: float
MJD of the observation.
mjd0: float
Reference MJD corresponding to (ra0, dec0). Default is J2000.
Examples
--------
>>> # Zero values have constant ra, dec
>>> get_ra_dec(0.0, 0.0, 0.0, 0.0, 0.0)
(np.float64(0.0), np.float64(0.0))
>>> mjd0 = 51544.5
>>> # 1 mas/yr for 1 year
>>> get_ra_dec(0.0, 0.0, 1.0, 1.0, mjd0 + 365.25, mjd0)
(np.float64(2.777777777455592e-07), np.float64(2.777777777455592e-07))
>>> # 1 mas/yr for 1 year near the pole
>>> get_ra_dec(0.0, 85.0, 1.0, 1.0, mjd0 + 365.25, mjd0)
(np.float64(3.1871427450005572e-06), np.float64(85.00000027777779))
"""
pm_ra_cosdec = np.where(np.isfinite(pm_ra_cosdec), pm_ra_cosdec, 0.0)
pm_dec = np.where(np.isfinite(pm_dec), pm_dec, 0.0)
t0 = Time(f"{mjd0}", format="mjd")
t1 = Time(f"{mjd}", format="mjd")
c0 = SkyCoord(
ra=ra0 * u.deg,
dec=dec0 * u.deg,
pm_ra_cosdec=pm_ra_cosdec * u.mas / u.yr,
pm_dec=pm_dec * u.mas / u.yr,
obstime=t0,
)
c1 = c0.apply_space_motion(new_obstime=t1)
return c1.ra.value, c1.dec.value
[docs]
def convert_flux(S1, E1_min=2, E1_max=10, E2_min=2, E2_max=7, Gamma=1.9):
"""
Convert flux S from bandpass E1 to bandpass E2.
Assumes a power-law spectrum with photon index Gamma.
Parameters
----------
S1: float
Input flux.
E1_min: float
Minimum energy in the input band.
E1_max: float
Maximum energy in the input band.
E2_min: float
Minimum energy in the output band.
E2_max: float
Maximum energy in the output band.
Gamma: float
Power-law photon index.
Returns
-------
S2: float
Converted flux in the output band.
Examples
--------
>>> # default band conversion
>>> convert_flux(1.0)
np.float64(0.7643018524251657)
>>> # modify Gamma
>>> convert_flux(1.0, Gamma=1.8)
np.float64(0.7498364916445219)
>>> # modify the maximum energy of the input band
>>> convert_flux(1.0, E1_max=8.0)
np.float64(0.8975323343244697)
>>> # Gamma=2.0 returns a non-finite value
>>> convert_flux(1.0, Gamma=2.0)
masked
"""
idx = 2 - Gamma
return S1 * np.ma.true_divide(E2_max**idx - E2_min**idx, E1_max**idx - E1_min**idx)
[docs]
def luminosity_to_flux(wavlen, luminosity_nu, redshift, distance_in_cm, use_igm=True):
"""
Convert luminosity_nu (in erg/s/Hz) to flux in uJy. Default distance is 10pc.
Parameters
----------
wavlen: float
Rest-frame wavelength in angstroms.
luminosity_nu: float
Rest-frame monochromatic luminosity in erg/s/Hz.
redshift: float
Redshift of the source.
distance_in_cm: float
Luminosity distance in cm.
use_igm: bool
Apply reddening by the intergalactic medium?
Returns
-------
flux: float
Flux in uJy at the given redshift.
Examples
--------
>>> from astropy.cosmology import FlatLambdaCDM
>>> import astropy.units as u
>>> cosmo = FlatLambdaCDM(H0=70.0, Om0=0.30)
>>> luminosity_to_flux(1.0, 1e32, 1.0, cosmo.luminosity_distance(1.0).cgs.value, True)
(0.00020000000000000004, np.float64(7.868437162608212e-16))
>>> luminosity_to_flux(1.0, 1e32, 1.0, cosmo.luminosity_distance(1.0).cgs.value, False)
(0.00020000000000000004, np.float64(1.2770363991236881e-15))
>>> # With 0 redshift the distance must be 10pc
>>> luminosity_to_flux(1.0, 1e32, 0.0, 0.0)
Traceback (most recent call last):
...
ValueError: For z=0, distance must correspond to 10pc.
>>> luminosity_to_flux(1.0, 1e32, 0.0, 10 * u.pc.to(u.cm), False)
(0.00010000000000000002, np.float64(278.78431938176107))
"""
if redshift == 0 and not np.isclose(distance_in_cm, (10 * u.pc).to(u.cm).value):
raise ValueError("For z=0, distance must correspond to 10pc.")
# Wavlen in angstrom and to observed frame
wavlen_observed = wavlen * (1 + redshift)
# To uJy in observed frame
log_flux = (
np.log10(luminosity_nu)
+ np.log10(wavlen_observed)
+ np.log10((u.erg / u.s / u.cm**2 / u.Hz).to(u.uJy))
- np.log10(constants.c.to(u.angstrom / u.s).value)
- np.log10(4 * np.pi)
- 2 * np.log10(distance_in_cm)
)
# Add IGM
if redshift > 0.0 and use_igm:
t_igm = my_get_IGM_absorption(redshift, lambda_obs=wavlen_observed)
log_flux += np.log10(t_igm)
# NOTE: everything is now in EGG units. Wavlen in um, flux in uJy
return wavlen_observed * (u.angstrom.to(u.um)), 10**log_flux
[docs]
def get_log_y_lo_hi(y, dy, null=99):
"""
Return logarithmic lower and upper limits assuming linear errors.
Examples
--------
>>> # Zero dy returns error
>>> get_log_y_lo_hi(0.0, 0.0)
(masked, masked, masked)
>>> # Test 10% relative error
>>> y0, y1, y2 = get_log_y_lo_hi(np.array([1.0]), np.array([0.10]))
>>> (y0.data, y1.data, y2.data)
(array([0.]), array([0.04575749]), array([0.04139269]))
"""
return (
np.ma.log10(y),
np.ma.log10(np.ma.true_divide(y, y - dy)),
np.ma.log10(np.ma.true_divide(y + dy, y)),
)
[docs]
def distance_modulus_to_parallax(mu):
"""
Convert distance module in mag to a parallax in mas.
Examples
--------
>>> distance_modulus_to_parallax(0.0)
np.float64(100.0)
>>> distance_modulus_to_parallax(1.0)
np.float64(63.09573444801932)
>>> distance_modulus_to_parallax(2.0)
np.float64(39.81071705534973)
"""
# NOTE: solved from mu \equiv 5 * log10(d) - 5
d = 10 ** (1 + mu / 5) * u.pc
return ((1 * u.au / d).si * u.rad).to(u.mas).value
[docs]
def _get_ratio_estimated_true(value_estimated: float, value_true: float) -> float:
"""
Calculate ratio between estimated value and true value.
Parameters
----------
value_estimated: float
Estimated value.
value_true: float
True value.
Returns
-------
ratio: float
The ratio defined as (value_estimated - value_true) / value_true.
Examples
--------
>>> _get_ratio_estimated_true(1.0, 1.0)
np.float64(0.0)
>>> _get_ratio_estimated_true(2.0, 1.0)
np.float64(1.0)
>>> _get_ratio_estimated_true(99.0, 0.0)
masked
"""
return np.ma.true_divide(np.abs(value_estimated - value_true), value_true)
[docs]
def get_sigma_nmad(value_estimated, value_true):
"""
Calculate sigma_NMAD from the given set of estimated / true values.
Reference is Hoaglin+ 1983. See also Sec. 4.1 of
https://iopscience.iop.org/article/10.1088/0004-637X/690/2/1236/meta
Parameters
----------
value_estimated: float
Estimated value.
value_true: float
True value.
Examples
--------
>>> get_sigma_nmad(1.0, 1.00)
np.float64(0.0)
>>> get_sigma_nmad(1.0, 0.10)
np.float64(13.32)
>>> get_sigma_nmad(1.0, 0.01)
np.float64(146.52)
"""
return 1.48 * np.median(_get_ratio_estimated_true(value_estimated, value_true))
[docs]
def get_fraction_catastrophic_error(value_estimated, value_true, limit=0.15):
"""
Calculate catastrophic error fraction from the set of estimated / true values.
Parameters
----------
value_estimated: float
Estimated value.
value_true: float
True value.
Examples
--------
>>> get_fraction_catastrophic_error(1.0, 1.0)
Traceback (most recent call last):
...
AttributeError: 'float' object has no attribute 'size'
>>> a = np.array([1, 2, 3])
>>> b = np.array([1, 1, 1])
>>> get_fraction_catastrophic_error(a, b)
np.float64(0.6666666666666666)
"""
n_total = value_estimated.size
is_catastrophic = _get_ratio_estimated_true(value_estimated, value_true) > limit
return is_catastrophic.sum() / n_total
[docs]
def get_log_lambda_SAR(i, N, m, z, t, seed):
"""
Parallelize Zou+2024 get_log_lambda_SAR.
Parameters
----------
i: int
UID of the source.
N: int
Total size of the catalog.
m: float
Host galaxy stellar mass.
z: float
Host galaxy redshift.
t: str
Host galaxy type.
seed: int
Random number seed.
Examples
--------
>>> get_log_lambda_SAR(0, 1, 9.5, 1.0, "star-forming", 222)
array(31.38098838)
"""
# NOTE: turns out that calling this function in parallel is probably not
# thread-safe. Enforce a fixed random number seed here
np.random.seed(seed)
ret = zou2024.get_log_lambda_SAR(m, z, t, add_ctk=True, log_mstar_lim=(9.50, 12.00), z_lim=(0.00, 4.00))
return ret
[docs]
def get_galaxy_ab(reff, ratio):
"""
Return the geometric mean-corrected galaxy semi-major and semi-major axes.
Starting with the definition of ellipticity = 1 - ratio, return 'a' and 'b'
defined through the geometric mean.
Parameters
----------
reff: float
Galaxy effective radius.
ratio: float
The ratio between the major and minor axes i.e. a/b.
Returns
-------
a, b: float
the 'a' and 'b' components defined as
(a, b) = (r_eff / sqrt(ratio), r_eff * sqrt(ratio))
Examples
--------
>>> get_galaxy_ab(1.0, 1.0)
(np.float64(1.0), np.float64(1.0))
>>> get_galaxy_ab(1.0, 0.5)
(np.float64(1.414213562373095), np.float64(0.7071067811865476))
"""
# ellipticity
# f = (a - b) / a
# f = 1 - b / a
# ratio = 1 -> perfectly round -> ellipticity = 0.0
# ratio = 0 -> perfectly round -> ellipticity = 1.0 (b = 0.0)
# ratio == b / a
# b = ratio * a
# ellip = 1 - ratio
a = np.ma.true_divide(reff, np.ma.sqrt(ratio))
b = reff * np.ma.sqrt(ratio)
return a, b
[docs]
def create_directory(filename: str) -> str:
"""
Create a directory corresponding to the filename.
If filename ends with '/', then filename is interpreted as the name of a
directory to be created. Otherwise the directory containing 'filename' is
created.
Parameters
----------
filename: str
Filename to create.
Returns
-------
dirname: str
Path to the directory that was created.
Raises
------
TypeError
If invalid type is given.
Examples
--------
>>> import os
>>> dirname = os.path.join("data", "tests", "test_util")
>>> if os.path.exists(dirname):
... os.rmdir(dirname)
>>> create_directory(f"{dirname}")
'data/tests'
>>> create_directory(f"{dirname}/")
'data/tests/test_util'
>>> create_directory(f"{dirname}/foo.bar")
'data/tests/test_util'
>>> create_directory(f"{dirname}/foo.bar/")
'data/tests/test_util/foo.bar'
>>> create_directory(f"{dirname}/bar.baz/test.dat")
'data/tests/test_util/bar.baz'
>>> create_directory(None)
Traceback (most recent call last):
...
TypeError: stat: path should be string, bytes, os.PathLike or integer, not NoneType
"""
dirname = filename
if not os.path.isdir(dirname):
dirname = os.path.dirname(dirname)
if not os.path.exists(dirname):
logger.info(f"Creating directory {dirname}")
os.makedirs(dirname)
if dirname.endswith("/"):
return dirname[:-1]
return dirname
[docs]
def get_mjd_vec():
"""
Return default MJD vector spanning ten-years with a delta of one day.
This is a simple convenience function to record the MJD vector in a single
function instead of a global variable.
Examples
--------
>>> get_mjd_vec()
array([ 0, 1, 2, ..., 3650, 3651, 3652], shape=(3653,))
"""
return np.arange(0, 3653, 1)
[docs]
def get_stellar_mass_completeness_cosmos2020(type: str, redshift: float) -> float:
"""
Return the 70% mass completeness limit according to COSMOS2020 SMF.
See https://arxiv.org/pdf/2212.02512 Eqs. 3, 4 and 5.
Returns
-------
stellar_mass_completeness: float or ArrayLike
70% stellar mass completeness limit in Msun
Examples
--------
>>> get_stellar_mass_completeness_cosmos2020("Total", 0.0)
46000000.0
>>> get_stellar_mass_completeness_cosmos2020("Total", 1.0)
248600000.0
>>> get_stellar_mass_completeness_cosmos2020("Star-forming", 1.0)
231000000.0
>>> get_stellar_mass_completeness_cosmos2020("non-existing type", 1.0)
Traceback (most recent call last):
...
KeyError: 'non-existing type'
"""
factors = {
"Total": (-3.23e7, 7.83e7),
"Star-forming": (-5.77e7, 8.66e7),
"Quiescent": (-3.79e7, 2.98e8),
}
f1, f2 = factors[type]
return f1 * (1 + redshift) + f2 * (1 + redshift) ** 2
[docs]
def read_fits(filename: str, *args, **kwargs):
"""
Read a FITS filename with supressed error messages.
"""
import warnings
from astropy.units import UnitsWarning
logger.info(f"Reading {filename} {args=} {kwargs=}")
ret = None
with warnings.catch_warnings():
warnings.simplefilter("ignore", UnitsWarning)
ret = fitsio.read(filename, *args, **kwargs)
return ret
[docs]
def read_table(filename: str, *args, **kwargs):
"""
Read an astropy table but do so silently.
"""
import warnings
from astropy.table import Table
from astropy.units import UnitsWarning
logger.info(f"Reading {filename} {args=} {kwargs=}")
with warnings.catch_warnings():
warnings.simplefilter("ignore", UnitsWarning)
table = Table.read(filename, *args, **kwargs)
return table
[docs]
def my_trapezoid(*args, **kwargs):
try:
return np.trapezoid(*args, **kwargs)
except AttributeError:
return np.trapz(*args, **kwargs)