Source code for gammapy.utils.random.utils

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Helper functions to work with distributions."""
import numbers
import numpy as np
import scipy.integrate
from astropy.coordinates import Angle
from astropy.time import TimeDelta

__all__ = [
    "draw",
    "get_random_state",
    "normalize",
    "pdf",
    "sample_powerlaw",
    "sample_sphere",
    "sample_sphere_distance",
    "sample_times",
]



[docs]
def normalize(func, x_min, x_max):
    """Normalize a 1D function over a given range."""

    def f(x):
        return func(x) / scipy.integrate.quad(func, x_min, x_max)[0]

    return f




[docs]
def pdf(func):
    """One-dimensional PDF of a given radial surface density."""

    def f(x):
        return x * func(x)

    return f




[docs]
def draw(low, high, size, dist, random_state="random-seed", *args, **kwargs):
    """Allow drawing of random numbers from any distribution."""
    from .inverse_cdf import InverseCDFSampler

    n = 1000
    x = np.linspace(low, high, n)

    pdf = dist(x)
    sampler = InverseCDFSampler(pdf=pdf, random_state=random_state)

    idx = sampler.sample(size)
    x_sampled = np.interp(idx, np.arange(n), x)
    return np.squeeze(x_sampled)




[docs]
def get_random_state(init):
    """Get a `numpy.random.RandomState` instance.

    The purpose of this utility function is to have a flexible way
    to initialise a `~numpy.random.RandomState` instance,
    a.k.a. a random number generator (``rng``).

    See :ref:`dev_random` for usage examples and further information.

    Parameters
    ----------
    init : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
        Available options to initialise the RandomState object:

        * ``int`` -- new RandomState instance seeded with this integer
          (calls `~numpy.random.RandomState` with ``seed=init``)
        * ``'random-seed'`` -- new RandomState instance seeded in a random way
          (calls `~numpy.random.RandomState` with ``seed=None``)
        * ``'global-rng'``, return the RandomState singleton used by ``numpy.random``.
        * `~numpy.random.RandomState` -- do nothing, return the input.

    Returns
    -------
    random_state : `~numpy.random.RandomState`
        RandomState instance.
    """
    if isinstance(init, (numbers.Integral, np.integer)):
        return np.random.RandomState(init)
    elif init == "random-seed":
        return np.random.RandomState(None)
    elif init == "global-rng":
        return np.random.mtrand._rand
    elif isinstance(init, np.random.RandomState):
        return init
    else:
        raise ValueError(
            "{} cannot be used to seed a numpy.random.RandomState"
            " instance".format(init)
        )




[docs]
def sample_sphere(size, lon_range=None, lat_range=None, random_state="random-seed"):
    """Sample random points on the sphere.

    Reference: http://mathworld.wolfram.com/SpherePointPicking.html

    Parameters
    ----------
    size : int
        Number of samples to generate.
    lon_range : `~astropy.coordinates.Angle`, optional
        Longitude range (min, max).
    lat_range : `~astropy.coordinates.Angle`, optional
        Latitude range (min, max).
    random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional
        Defines random number generator initialisation.
        Passed to `~gammapy.utils.random.get_random_state`.
        Default is "random-seed".

    Returns
    -------
    lon, lat: `~astropy.coordinates.Angle`
        Longitude and latitude coordinate arrays.
    """
    random_state = get_random_state(random_state)

    if lon_range is None:
        lon_range = Angle([0.0, 360.0], "deg")

    if lat_range is None:
        lat_range = Angle([-90.0, 90.0], "deg")

    # Sample random longitude
    u = random_state.uniform(size=size)
    lon = lon_range[0] + (lon_range[1] - lon_range[0]) * u

    # Sample random latitude
    v = random_state.uniform(size=size)
    z_range = np.sin(lat_range)
    z = z_range[0] + (z_range[1] - z_range[0]) * v
    lat = np.arcsin(z)

    return lon, lat




[docs]
def sample_sphere_distance(
    distance_min=0, distance_max=1, size=None, random_state="random-seed"
):
    """Sample random distances if the 3-dim space density is constant.

    This function uses inverse transform sampling
    (`Wikipedia <http://en.wikipedia.org/wiki/Inverse_transform_sampling>`__)
    to generate random distances for an observer located in a 3-dim
    space with constant source density in the range ``(distance_min, distance_max)``.

    Parameters
    ----------
    distance_min : float, optional
        Minimum distance. Default is 0.
    distance_max : float, optional
        Maximum distance. Default is 1.
    size : int or tuple of ints, optional
        Output shape. Default is one sample. Passed to `numpy.random.uniform`.
    random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional
        Defines random number generator initialisation.
        Passed to `~gammapy.utils.random.get_random_state`.
        Default is "random-seed".

    Returns
    -------
    distance : `~numpy.ndarray`
        Array of samples.
    """
    random_state = get_random_state(random_state)

    # Since the differential distribution is dP / dr ~ r ^ 2,
    # we have a cumulative distribution
    #     P(r) = a * r ^ 3 + b
    # with P(r_min) = 0 and P(r_max) = 1 implying
    #     a = 1 / (r_max ^ 3 - r_min ^ 3)
    #     b = -a * r_min ** 3

    a = 1.0 / (distance_max**3 - distance_min**3)
    b = -a * distance_min**3

    # Now for inverse transform sampling we need to use the inverse of
    #     u = a * r ^ 3 + b
    # which is
    #     r = [(u - b)/ a] ^ (1 / 3)
    u = random_state.uniform(size=size)
    distance = ((u - b) / a) ** (1.0 / 3)

    return distance




[docs]
def sample_powerlaw(x_min, x_max, gamma, size=None, random_state="random-seed"):
    """Sample random values from a power law distribution.

    f(x) = x ** (-gamma) in the range x_min to x_max

    It is assumed that *gamma* is the **differential** spectral index.

    Reference: http://mathworld.wolfram.com/RandomNumber.html

    Parameters
    ----------
    x_min : float
        Minimum x range.
    x_max : float
        Maximum x range.
    gamma : float
        Power law index.
    size : int, optional
        Number of samples to generate. Default is None.
    random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional
        Defines random number generator initialisation.
        Passed to `~gammapy.utils.random.get_random_state`.
        Default is "random-seed".

    Returns
    -------
    x : `~numpy.ndarray`
        Array of samples from the distribution.
    """
    random_state = get_random_state(random_state)

    size = int(size)

    exp = -gamma
    base = random_state.uniform(x_min**exp, x_max**exp, size)
    x = base ** (1 / exp)

    return x




[docs]
def sample_times(
    size,
    rate,
    dead_time=TimeDelta(0, format="sec"),
    return_diff=False,
    random_state="random-seed",
):
    """Make random times assuming a Poisson process.

    This function can be used to test event time series,
    to have a comparison what completely random data looks like.

    Can be used in two ways (in either case the return type is `~astropy.time.TimeDelta`):

    * ``return_delta=False`` - Return absolute times, relative to zero (default)
    * ``return_delta=True`` - Return time differences between consecutive events.

    Parameters
    ----------
    size : int
        Number of samples.
    rate : `~astropy.units.Quantity`
        Event rate.
    dead_time : `~astropy.units.Quantity` or `~astropy.time.TimeDelta`, optional
        Dead time after event.
    return_diff : bool
        Return time difference between events. Default False, return absolute times.
    random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}, optional
        Defines random number generator initialisation.
        Passed to `~gammapy.utils.random.get_random_state`.
        Default is "random-seed".

    Returns
    -------
    time : `~astropy.time.TimeDelta`
        Time differences (second) after time zero.

    Examples
    --------
    Example how to simulate 100 events at a rate of 10 Hz.
    As expected the last event occurs after about 10 seconds.

    >>> from astropy.units import Quantity
    >>> from gammapy.utils.random import sample_times
    >>> rate = Quantity(10, 'Hz')
    >>> times = sample_times(size=100, rate=rate, random_state=0)
    >>> times[-1]
    <TimeDelta object: scale='None' format='sec' value=9.186484131475074>
    """
    random_state = get_random_state(random_state)

    dead_time = TimeDelta(dead_time)
    scale = (1 / rate).to("s").value
    time_delta = random_state.exponential(scale=scale, size=size)
    time_delta += dead_time.to("s").value

    if return_diff:
        return TimeDelta(time_delta, format="sec")
    else:
        time = time_delta.cumsum()
        return TimeDelta(time, format="sec")