increase util.py doctest coverage to 100% (e90ff718) · Commits · Akke Esa Tapio Viitanen / lsst_inaf_agile

src/lsst_inaf_agile/util.py

+441 −90

Original line number	Diff line number	Diff line
		@@ -25,45 +25,222 @@ logger = logging.getLogger(__name__)


		def flux_to_mag(flux):
		"""Convert uJy flux to AB magnitude."""
		"""
		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()


		def mag_to_flux(mag):
		"""Convert AB magnitude to uJy flux."""
		"""
		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)


		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(mag))
		flux = np.sum(mag_to_flux(np.array(mag)))
		return flux_to_mag(flux)


		def get_volume(zmin, zmax, area_deg2, H0=70.0, Om0=0.30, Tcmb0=2.73):
		"""Return the comoving volume in Mpc for a redshift shell."""
		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
		Om0: float
		Present day dimensionless matter density parameter.
		Tcmb0: float
		Present day CMB temperature in Kelvin.

		Examples
		--------
		The following values only illustrate the usage. The accuracy of the volume
		estimation is set by AstroPy 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


		def get_chisq_nu(y, y_model, sigma):
		"""Return the chisq value of the measurement y."""
		diff = (y - y_model) / sigma
		return np.sum(diff**2) / (y.size - 1)
		"""
		Return the chisq value of the measurement y.

		Input arguments are converted to numpy arrays before calculation.

		def get_key_function(bins, x, values=None, nmin=30, args, *kwargs):
		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)


		def get_key_function(bins, x, values, 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)
		@@ -90,14 +267,65 @@ def get_key_function(bins, x, values=None, nmin=30, args, *kwargs):
		return centers, dbins / 2, function, dfunction


		def egg_band_to_index(egg, band: str) -> int:
		"""Convert EGG band name to an index."""
		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'.

		Examples
		--------
		>>> # Generate a mock EGG catalog
		>>> 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
		"""
		bands = [b.strip() for b in egg["BANDS"][0]]
		return bands.index(band)


		def get_ra_dec(ra0, dec0, pm_ra_cosdec, pm_dec, mjd, mjd0=51544.5):
		"""Get ra, dec for current epoch."""
		"""
		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)

		@@ -120,74 +348,81 @@ 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.
		"""
		idx = 2 - Gamma
		return S1 * np.true_divide(E2_maxidx - E2_minidx, E1_maxidx - E1_minidx)

		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.

		def get_log_L_2_keV(log_LX_2_10, Gamma=1.9, wavelength=6.2):
		"""
		Return monochromatic X-ray luminosity at lambda = wavelength in erg/s Hz^-1.
		Returns
		-------
		S2: float
		Converted flux in the output band.

		To be used for the alpha_ox Lx = restframe 2-10 kev luminosity.
		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
		"""
		Lx = 10**log_LX_2_10
		K = (Lx / (6.2 (Gamma - 2) - 1.24 (Gamma - 2))) * (Gamma - 2) # 6.2, 1.24 = 2kev, 10kev in A°
		return np.log10((K * wavelength ** (Gamma - 1)) / 2.998e18)
		idx = 2 - Gamma
		return S1 * np.ma.true_divide(E2_maxidx - E2_minidx, E1_maxidx - E1_minidx)


		def get_log_L_2500(log_L_2_keV, alpha=0.952, beta=2.138, scatter=True):
		def luminosity_to_flux(wavlen, luminosity, redshift, distance_in_cm, use_igm=True):
		"""
		Return the 2500 ang° monochromatic luminosity (in erg/s).

		It uses Lusso+10 eq. 5 (inverted) Lx = alpha L_opt - beta.
		Convert luminosity (in erg/s/ang) to flux in uJy. Default distance is 10pc.

		Parameters
		----------
		wavlen: float
		Rest-frame wavelength in angstroms.
		luminosity: float
		Rest-frame luminosity in erg/s/angstrom.
		redshift: float
		Redshift of the source.
		distance_in_cm: float
		Luminosity distance in cm.
		use_igm: bool
		Apply reddening by the intergalactic medium?
		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))
		"""
		log_L_2500 = (log_L_2_keV + beta) / alpha
		assert np.allclose(alpha * log_L_2500 - beta, log_L_2_keV)

		# TODO: implement realistic scatter
		if scatter:
		log_L_2500 += np.random.normal(loc=0, scale=0.4, size=log_L_2500.size)

		return log_L_2500


		def get_E_BV(
		type2=False,
		alpha_1=7.93483055,
		n_1=2.97565676,
		alpha_2=11.6133635,
		n_2=1.42972,
		mu_type_2=0.3,
		):
		"""Return E(B-V) using the functional form from Hopkins+2004."""
		type_1_ebv = (np.linspace(0, 1, 101),)
		type_2_ebv = (np.linspace(0, 3, 301),)

		def sample_ebv(N_AGN, probability_distribution, ebv_range, *args):
		"""Sample the E(B-V) distribution."""
		cumulative = np.cumsum(probability_distribution(ebv_range, *args))
		cumulative /= np.max(cumulative)
		return np.interp(np.random.rand(N_AGN), cumulative, ebv_range)

		def hopkins04(x, alpha, n):
		"""Return p(E_BV)."""
		y = 1 / (1 + (x * alpha) ** n)
		return y / np.trapezoid(y, x)

		ebv = None
		if type2:
		ebv = sample_ebv(1, hopkins04, type_2_ebv, alpha_2, n_2) + mu_type_2
		else:
		ebv = sample_ebv(1, hopkins04, type_1_ebv, alpha_1, n_1)
		return np.squeeze(ebv)


		def luminosity_to_flux(wavlen, luminosity, redshift, distance_in_cm, use_igm=True):
		"""Convert luminosity (in erg/s/ang) to flux in uJy. Default distance is 10pc."""
		if redshift == 0:
		assert np.isclose(distance_in_cm, (10 * u.pc).to(u.cm).value)
		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)
		@@ -212,22 +447,42 @@ def luminosity_to_flux(wavlen, luminosity, redshift, distance_in_cm, use_igm=Tru


		def get_log_y_lo_hi(y, dy, null=99):
		"""Return logarithmic lower and upper limits assuming linear errors."""
		y0 = np.ma.log10(y)
		y1 = np.ma.log10(y / (y - dy))
		y2 = np.ma.log10((y + dy) / y)

		select = y - dy <= 0.0
		y1[select] = null
		"""
		Return logarithmic lower and upper limits assuming linear errors.

		return y0, y1, y2
		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)),
		)


		def distance_modulus_to_parallax(mu):
		"""Convert distance module to a parallax."""
		"""
		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
		B = 1 * u.au
		return ((B / d).si * u.rad).to(u.mas).value
		return ((1 * u.au / d).si * u.rad).to(u.mas).value


		def get_star_binary_fbin(star, binary, fbin=0.40, nrepeat=4, seed=1206):
		@@ -247,6 +502,20 @@ def get_star_binary_fbin(star, binary, fbin=0.40, nrepeat=4, seed=1206):
		number of binary systems present in the same regions.
		"

		Examples
		--------
		>>> from astropy.table import Table
		>>> star = Table({"a": [0.0, 1.0, 2.0] * 1000})
		>>> binary = Table({"a": [0.0, 1.0, 2.0] * 100, "b": [0.0, 1.0, 2.0] * 100})
		>>> star2, binary2 = get_star_binary_fbin(star, binary)
		>>> len(star2), len(binary2)
		(1785, 1200)
		>>> # Small catalog will warn about insufficient statistics but succeeds
		>>> star = Table({"a": [0.0, 1.0, 2.0] * 5})
		>>> binary = Table({"a": [0.0, 1.0, 2.0] * 2, "b": [0.0, 1.0, 2.0] * 2})
		>>> star2, binary2 = get_star_binary_fbin(star, binary)
		>>> len(star2), len(binary2)
		(10, 24)
		"""
		# Set the seed
		np.random.seed(seed)
		@@ -258,13 +527,12 @@ def get_star_binary_fbin(star, binary, fbin=0.40, nrepeat=4, seed=1206):
		binary2 = np.repeat(binary, nrepeat)

		# Copy over ra/dec/etc from the REMAINING stellar catalog

		# NOTE: in some small catalog cases the number of remaining stars is not
		# enough to sample for the binary catalog. In these cases, set replace=True
		is_not_enough = len(binary2) > (~is_star).sum()
		if is_not_enough:
		logger.warning(
		"Small stellar catalog. Can not sample binary stars sufficiently.Will use replace=True"
		"Small stellar catalog. " "Can not sample binary stars sufficiently. " "Will use replace=True"
		)
		star2 = np.random.choice(star[~is_star], size=len(binary2), replace=is_not_enough)

		@@ -280,6 +548,13 @@ def _get_ratio_estimated_true(value_estimated: float, value_true: float) -> floa
		Calculate ratio between estimated value and true value.

		The "ratio" is defined as (y_est - y_true) / y_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)
		"""
		return np.ma.true_divide(np.abs(value_estimated - value_true), value_true)

		@@ -289,14 +564,35 @@ 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

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


		def get_fraction_catastrophic_error(value_estimated, value_true, limit=0.15):
		"""Calculate catastrophic error fraction from the set of estimated / true values."""
		"""
		Calculate catastrophic error fraction from the set of estimated / true values.

		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
		@@ -320,6 +616,11 @@ def get_log_lambda_SAR(i, N, m, z, t, seed):
		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
		@@ -343,6 +644,12 @@ def get_galaxy_ab(reff, ratio):
		b: float
		the 'b' component: 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
		@@ -388,7 +695,18 @@ def create_directory(filename: str) -> None:


		def get_mjd_vec():
		"""Return default MJD vector spanning ten-years with a delta of one day."""
		"""
		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)


		@@ -403,6 +721,18 @@ def get_stellar_mass_completeness_cosmos2020(type: str, redshift: float) -> floa
		stellar_mass_completeness: float or array_like
		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),
		@@ -415,7 +745,15 @@ def get_stellar_mass_completeness_cosmos2020(type: str, redshift: float) -> floa

		def read_fits(filename, args, *kwargs):
		"""
		Read a FITS filename with supressed error messages
		Read a FITS filename with supressed error messages.

		Examples
		--------
		>>> import fitsio
		>>> fitsio.write("my_fits_file.fits", {"a": np.array([0, 1, 2])}, clobber=True)
		>>> read_fits("my_fits_file.fits")
		array([(0,), (1,), (2,)], dtype=[('a', '>i8')])
		>>> os.remove("my_fits_file.fits")
		"""
		import warnings

		@@ -430,7 +768,20 @@ def read_fits(filename, args, *kwargs):


		def read_table(filename):
		"""Read an astropy table but do so silently."""
		"""
		Read an astropy table but do so silently.

		Examples
		--------
		>>> import fitsio
		>>> fitsio.write("my_fits_file.fits", {"a": np.array([0, 1, 2])}, clobber=True)
		>>> read_table("my_fits_file.fits")["a"]
		<Column name='a' dtype='int64' length=3>
		0
		1
		2
		>>> os.remove("my_fits_file.fits")
		"""
		import warnings

		from astropy.table import Table