lsst_inaf_agile.util

General utilities.

Attributes

logger

Functions

flux_to_mag(flux)	Convert uJy flux to AB magnitude.
mag_to_flux(mag)	Convert AB magnitude to uJy flux.
mag_sum(mag)	Sum together magnitudes 'mag', return the combined magnitude.
get_volume(zmin, zmax[, area_deg2, H0, Om0, Tcmb0])	Return the flat LambdaCDM comoving volume in Mpc for a redshift shell.
get_chisq_nu(y, y_model, sigma)	Return the chisq value of the measurement y.
get_key_function(bins, x[, values, nmin])	Return the "key" function i.e. number of objects at interval in 'key'.
egg_band_to_index(→ int)	Convert EGG band name to an index.
get_ra_dec(ra0, dec0, pm_ra_cosdec, pm_dec, mjd[, mjd0])	Get ra, dec for current epoch modified by the proper motion.
convert_flux(S1[, E1_min, E1_max, E2_min, E2_max, Gamma])	Convert flux S from bandpass E1 to bandpass E2.
luminosity_to_flux(wavlen, luminosity_nu, redshift, ...)	Convert luminosity_nu (in erg/s/Hz) to flux in uJy. Default distance is 10pc.
get_log_y_lo_hi(y, dy[, null])	Return logarithmic lower and upper limits assuming linear errors.
distance_modulus_to_parallax(mu)	Convert distance module in mag to a parallax in mas.
_get_ratio_estimated_true(→ float)	Calculate ratio between estimated value and true value.
get_sigma_nmad(value_estimated, value_true)	Calculate sigma_NMAD from the given set of estimated / true values.
get_fraction_catastrophic_error(value_estimated, ...)	Calculate catastrophic error fraction from the set of estimated / true values.
get_log_lambda_SAR(i, N, m, z, t, seed)	Parallelize Zou+2024 get_log_lambda_SAR.
get_galaxy_ab(reff, ratio)	Return the geometric mean-corrected galaxy semi-major and semi-major axes.
create_directory(→ str)	Create a directory corresponding to the filename.
get_mjd_vec()	Return default MJD vector spanning ten-years with a delta of one day.
get_stellar_mass_completeness_cosmos2020(→ float)	Return the 70% mass completeness limit according to COSMOS2020 SMF.
read_fits(filename, *args, **kwargs)	Read a FITS filename with supressed error messages.
read_table(filename)	Read an astropy table but do so silently.

Module Contents

logger

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)

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'

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)

get_volume(zmin: float, zmax: float, area_deg2: float = 41252.96124941928, H0: float = 70.0, Om0: float = 0.3, 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 – The comoving volume segment in Mpc3 corresponding to the input arguments.

Return type:

float

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)

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 – Estimated reduced chisq value. Degrees of freedom are assumed to be len(y)-1.

Return type:

float

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)

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 – Additional arguments forwarded to the function ‘get_volume’.

• kwargs – Additional arguments forwarded to the function ‘get_volume’.

Returns:

Bin centers (x), Bin width (dx), Counts (y), Delta counts (dy). Error on the counts is assumed to be Poissonian i.e. sqrt(Ncounts).

Return type:

x, dx, y, dy

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)

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 – Index of the band in the EGG catalog.

Return type:

int

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

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))

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 – Converted flux in the output band.

Return type:

float

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

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 – Flux in uJy at the given redshift.

Return type:

float

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))

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]))

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)

_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 – The ratio defined as (value_estimated - value_true) / value_true.

Return type:

float

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

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)

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)

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)

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 – the ‘a’ and ‘b’ components defined as (a, b) = (r_eff / sqrt(ratio), r_eff * sqrt(ratio))

Return type:

float

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))

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 – Path to the directory that was created.

Return type:

str

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

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,))

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 – 70% stellar mass completeness limit in Msun

Return type:

float or ArrayLike

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'

read_fits(filename, *args, **kwargs)

Read a FITS filename with supressed error messages.

read_table(filename)

Read an astropy table but do so silently.