Source code for halotools.empirical_models.model_helpers

This module contains general purpose helper functions
used by many of the Halotools models.

import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline as spline
from scipy.special import gammaincc, gamma, expi

from ..utils.array_utils import custom_len
from ..custom_exceptions import HalotoolsError

__all__ = ('solve_for_polynomial_coefficients', 'polynomial_from_table',
            'enforce_periodicity_of_box', 'custom_spline', 'create_composite_dtype',
            'bind_default_kwarg_mixin_safe', 'custom_incomplete_gamma')

__author__ = ['Andrew Hearin', 'Surhud More', 'Johannes Ulf Lange']

[docs] def solve_for_polynomial_coefficients(abscissa, ordinates): r""" Solves for coefficients of the unique, minimum-degree polynomial that passes through the input abscissa and attains values equal the input ordinates. Parameters ---------- abscissa : array Elements are the abscissa at which the desired values of the polynomial have been tabulated. ordinates : array Elements are the desired values of the polynomial when evaluated at the abscissa. Returns ------- polynomial_coefficients : array Elements are the coefficients determining the polynomial. Element i of polynomial_coefficients gives the degree i polynomial coefficient. Notes -------- Input arrays abscissa and ordinates can in principle be of any dimension Ndim, and there will be Ndim output coefficients. The input ordinates specify the desired values of the polynomial when evaluated at the Ndim inputs specified by the input abscissa. There exists a unique, order Ndim polynomial that returns the input ordinates when the polynomial is evaluated at the input abscissa. The coefficients of that unique polynomial are the output of the function. As an example, suppose that a model in which the quenched fraction is :math:`F_{q}(logM_{\mathrm{halo}} = 12) = 0.25` and :math:`F_{q}(logM_{\mathrm{halo}} = 15) = 0.9`. Then this function takes [12, 15] as the input abscissa, [0.25, 0.9] as the input ordinates, and returns the array :math:`[c_{0}, c_{1}]`. The unique polynomial linear in :math:`log_{10}M` that passes through the input ordinates and abscissa is given by :math:`F(logM) = c_{0} + c_{1}*log_{10}logM`. Examples -------- >>> abscissa = [0, 2] >>> ordinates = [0, 2] >>> coeff = solve_for_polynomial_coefficients(abscissa, ordinates) >>> assert np.allclose(coeff, (0, 1)) """ abscissa = np.atleast_1d(abscissa) ordinates = np.atleast_1d(ordinates) columns = np.ones(len(abscissa)) for i in np.arange(len(abscissa)-1): columns = np.append(columns, [abscissa**(i+1)]) quenching_model_matrix = columns.reshape( len(abscissa), len(abscissa)).transpose() polynomial_coefficients = np.linalg.solve( quenching_model_matrix, ordinates) return np.array(polynomial_coefficients)
[docs] def polynomial_from_table(table_abscissa, table_ordinates, input_abscissa): r""" Method to evaluate an input polynomial at the input_abscissa. The input polynomial is determined by `solve_for_polynomial_coefficients` from table_abscissa and table_ordinates. Parameters ---------- table_abscissa : array Elements are the abscissa determining the input polynomial. table_ordinates : array Elements are the desired values of the input polynomial when evaluated at table_abscissa input_abscissa : array Points at which to evaluate the input polynomial. Returns ------- output_ordinates : array Values of the input polynomial when evaluated at input_abscissa. Examples --------- >>> table_abscissa = [0, 1, 2, 3] >>> table_ordinates = [0, 2, 4, 6] >>> input_abscissa = 0.5 >>> result = polynomial_from_table(table_abscissa, table_ordinates, input_abscissa) >>> assert np.allclose(result, 1.0) """ input_abscissa = np.atleast_1d(input_abscissa) coefficient_array = solve_for_polynomial_coefficients( table_abscissa, table_ordinates) output_ordinates = np.zeros(custom_len(input_abscissa)) # Use coefficients to compute values of the inflection function polynomial for n, coeff in enumerate(coefficient_array): output_ordinates += coeff*input_abscissa**n return output_ordinates
[docs] def enforce_periodicity_of_box(coords, box_length, check_multiple_box_lengths=False, **kwargs): r""" Function used to apply periodic boundary conditions of the simulation, so that mock galaxies all lie in the range [0, Lbox]. Parameters ---------- coords : array_like float or ndarray containing a set of points with values ranging between [-box_length, 2*box_length] box_length : float the size of simulation box (currently hard-coded to be Mpc/h units) velocity : array_like, optional velocity in the same dimension as the input coords. check_multiple_box_lengths : bool, optional If True, an exception will be raised if the points span a range of more than 2Lbox. Default is False. Returns ------- periodic_coords : array_like array with values and shape equal to input coords, but with periodic boundary conditions enforced """ if check_multiple_box_lengths is True: xmin = np.min(coords) if xmin < -box_length: msg = ("\nThere is at least one input point with a coordinate less than -Lbox\n") raise HalotoolsError(msg) xmax = np.max(coords) if xmax > 2*box_length: msg = ("\nThere is at least one input point with a coordinate greater than 2*Lbox\n") raise HalotoolsError(msg) try: velocity = kwargs['velocity'] newcoords = coords % box_length return newcoords, velocity except: return coords % box_length
[docs] def custom_spline(table_abscissa, table_ordinates, **kwargs): r""" Convenience wrapper around `~scipy.interpolate.InterpolatedUnivariateSpline`, written specifically to handle the edge case of a spline table being built from a single point. Parameters ---------- table_abscissa : array_like abscissa values defining the interpolation table_ordinates : array_like ordinate values defining the interpolation k : int, optional Degree of the desired spline interpolation. Default is 1. Returns ------- output : object Function object to use to evaluate the interpolation of the input table_abscissa & table_ordinates Notes ----- Only differs from `~scipy.interpolate.UnivariateSpline` in two respects. First, the degree of the spline interpolation is automatically chosen to be the maximum allowable degree permitted by the number of abscissa points. Second, the behavior differs for the case where the input tables have only a single element. In this case, the default behavior of the scipy function is to raise an exception. The `custom_spline` instead returns a constant-valued function where the returned value is simply the scalar value of the input ordinates. """ if custom_len(table_abscissa) != custom_len(table_ordinates): len_abscissa = custom_len(table_abscissa) len_ordinates = custom_len(table_ordinates) raise HalotoolsError("table_abscissa and table_ordinates must have the same length \n" " len(table_abscissa) = %i and len(table_ordinates) = %i" % (len_abscissa, len_ordinates)) max_scipy_spline_degree = 5 if 'k' in kwargs: k = np.min([custom_len(table_abscissa)-1, kwargs['k'], max_scipy_spline_degree]) else: k = 1 if k < 0: raise HalotoolsError("Spline degree must be non-negative") elif k == 0: if custom_len(table_ordinates) != 1: raise HalotoolsError("In spline_degree=0 edge case, " "table_abscissa and table_abscissa must be 1-element arrays") return lambda x: np.zeros(custom_len(x)) + table_ordinates[0] else: spline_function = spline(table_abscissa, table_ordinates, k=k) return spline_function
def call_func_table(func_table, abscissa, func_indices): r""" Returns the output of an array of functions evaluated at a set of input points if the indices of required functions is known. Parameters ---------- func_table : array_like Length k array of function objects abscissa : array_like Length Npts array of points at which to evaluate the functions. func_indices : array_like Length Npts array providing the indices to use to choose which function operates on each abscissa element. Thus func_indices is an array of integers ranging between 0 and k-1. Returns ------- out : array_like Length Npts array giving the evaluation of the appropriate function on each abscissa element. """ func_table = np.atleast_1d(func_table) shape_error_msg = "Input ``func_table must be one-dimensional, but has shape = {0}" assert len(np.shape(func_table)) == 1, shape_error_msg.format(func_table.shape) abscissa = np.atleast_1d(abscissa) func_indices = np.atleast_1d(func_indices) func_argsort = func_indices.argsort() func_ranges = list(np.searchsorted(func_indices[func_argsort], list(range(len(func_table))))) func_ranges.append(None) out = np.zeros_like(abscissa) for f, start, end in zip(func_table, func_ranges[:-1], func_ranges[1:]): ix = func_argsort[start:end] out[ix] = f(abscissa[ix]) return out def bind_required_kwargs(required_kwargs, obj, **kwargs): r""" Method binds each element of ``required_kwargs`` to the input object ``obj``, or raises and exception for cases where a mandatory keyword argument was not passed to the ``obj`` constructor. Used throughout the package when a required keyword argument has no obvious default value. Parameters ---------- required_kwargs : list List of strings of the keyword arguments that are required when instantiating the input ``obj``. obj : object The object being instantiated. Notes ----- The `bind_required_kwargs` method assumes that each required keyword argument should be bound to ``obj`` as attribute with the same name as the keyword. """ for key in required_kwargs: if key in list(kwargs.keys()): setattr(obj, key, kwargs[key]) else: class_name = obj.__class__.__name__ msg = ( key + ' is a required keyword argument ' + 'to instantiate the '+class_name+' class' ) raise KeyError(msg)
[docs] def create_composite_dtype(dtype_list): r""" Find the union of the dtypes in the input list, and return a composite dtype after verifying consistency of typing of possibly repeated fields. Parameters ---------- dtype_list : list List of dtypes with possibly repeated field names. Returns -------- composite_dtype : dtype Numpy dtype object composed of the union of the input dtypes. Notes ----- Basically an awkward workaround to the fact that numpy dtype objects are not iterable. """ name_list = list(set([name for d in dtype_list for name in d.names])) composite_list = [] for name in name_list: for dt in dtype_list: if name in dt.names: tmp = np.dtype(composite_list) if name in tmp.names: if tmp[name].type == dt[name].type: pass else: msg = ("Inconsistent dtypes for name = ``%s``.\n" " dtype1 = %s\n dtype2 = %s\n" % (name, tmp[name].type, dt[name].type)) raise HalotoolsError(msg) else: composite_list.append((name, dt[name].type)) composite_dtype = np.dtype(composite_list) return composite_dtype
[docs] def bind_default_kwarg_mixin_safe(obj, keyword_argument, constructor_kwargs, default_value): r""" Function used to ensure that a keyword argument passed to the constructor of an orthogonal mix-in class is not already an attribute bound to self. If it is safe to bind the keyword_argument to the object, `bind_default_kwarg_mixin_safe` will do so. Parameters ---------- obj : class instance Instance of the class to which we want to bind the input ``keyword_argument``. keyword_argument : string name of the attribute that will be bound to the object if the action is deemed mix-in safe. constructor_kwargs : dict keyword argument dictionary passed to the constructor of the input ``obj``. default_value : object Whatever the default value for the attribute should be if ``keyword_argument`` does not appear in kwargs nor is it already bound to the ``obj``. """ if hasattr(obj, keyword_argument): if keyword_argument in constructor_kwargs: clname = obj.__class__.__name__ msg = ("Do not pass the ``%s`` keyword argument " "to the constructor of the %s class \nwhen using the %s class " "as an orthogonal mix-in" % (keyword_argument, clname, clname)) raise HalotoolsError(msg) else: pass else: if keyword_argument in constructor_kwargs: setattr(obj, keyword_argument, constructor_kwargs[keyword_argument]) else: setattr(obj, keyword_argument, default_value)
[docs] def custom_incomplete_gamma(a, x): r""" Incomplete gamma function. For the case covered by scipy, a > 0, scipy is called. Otherwise the gamma function recurrence relations are called, extending the scipy behavior. Parameters ----------- a : array_like x : array_like Returns -------- gamma : array_like Examples -------- >>> a, x = 1, np.linspace(1, 10, 100) >>> g = custom_incomplete_gamma(a, x) >>> a = 0 >>> g = custom_incomplete_gamma(a, x) >>> a = -1 >>> g = custom_incomplete_gamma(a, x) """ if isinstance(a, np.ndarray): if not isinstance(x, np.ndarray): x = np.repeat(x, len(a)) if len(a) != len(x): msg = ("The ``a`` and ``x`` arguments of the " "``custom_incomplete_gamma`` function must have the same" "length.\n") raise HalotoolsError(msg) result = np.zeros(len(a)) mask = (a < 0) if np.any(mask): result[mask] = ((custom_incomplete_gamma(a[mask]+1, x[mask]) - x[mask]**a[mask] * np.exp(-x[mask])) / a[mask]) mask = (a == 0) if np.any(mask): result[mask] = -expi(-x[mask]) mask = a > 0 if np.any(mask): result[mask] = gammaincc(a[mask], x[mask]) * gamma(a[mask]) return result else: if a < 0: return (custom_incomplete_gamma(a+1, x) - x**a * np.exp(-x))/a elif a == 0: return -expi(-x) else: return gammaincc(a, x) * gamma(a)
custom_incomplete_gamma.__author__ = ['Surhud More', 'Johannes Ulf Lange']