# Copyright (c) 2012-2023 by the GalSim developers team on GitHub
# https://github.com/GalSim-developers
#
# This file is part of GalSim: The modular galaxy image simulation toolkit.
# https://github.com/GalSim-developers/GalSim
#
# GalSim is free software: redistribution and use in source and binary forms,
# with or without modification, are permitted provided that the following
# conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice, this
# list of conditions, and the disclaimer given in the accompanying LICENSE
# file.
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions, and the disclaimer given in the documentation
# and/or other materials provided with the distribution.
#
__all__ = [ 'BaseDeviate', 'UniformDeviate', 'GaussianDeviate', 'PoissonDeviate', 'DistDeviate',
'BinomialDeviate', 'Chi2Deviate', 'GammaDeviate', 'WeibullDeviate', ]
import numpy as np
import os
from . import _galsim
from .errors import GalSimRangeError, GalSimValueError, GalSimIncompatibleValuesError
from .errors import galsim_warn
from ._utilities import isinteger, math_eval
from .table import LookupTable
from . import integ
[docs]class BaseDeviate:
"""Base class for all the various random deviates.
This holds the essential random number generator that all the other classes use.
All deviates take an initial ``seed`` argument that is used to seed the underlying random number
generator. It has three different kinds of behavior.
1. An integer value can be provided to explicitly seed the random number generator with a
particular value. This is useful to have deterministic behavior. If you seed with an
integer value, the subsequent series of "random" values will be the same each time you
run the program.
2. A seed of 0 or None means to pick some arbitrary value that will be different each time
you run the program. Currently, this tries to get a seed from /dev/urandom if possible.
If that doesn't work, then it creates a seed from the current time. You can also get this
behavior by omitting the seed argument entirely. (i.e. the default is None.)
3. Providing another BaseDeviate object as the seed will make the new BaseDeviate share the
same underlying random number generator as the other BaseDeviate. So you can make one
BaseDeviate (of any type), and seed it with a particular deterministic value. Then if you
pass that BaseDeviate to any other one you make, they will both be using the same RNG and
the series of "random" values will be deterministic.
**Usage**:
The only kind of random number you can obtain from a pure BaseDeviate (i.e. one that is
not actually one of the variosu subclasses) is a "raw" value. This is an unsigned 32-bit
integer that behind the scenes is used by all sub-classes to generate floating point values
with various distributions.
>>> rng = galsim.BaseDeviate(215324)
>>> rng.raw()
3559052779
Most other usage is effected by creating sub-classes corresponding to particular random
deviates with various distributions. E.g. `UniformDeviate` generates random values following
a uniform distribution between 0 and 1.
>>> rng = galsim.BaseDeviate(215324)
>>> ud = galsim.UniformDeviate(rng)
>>> ud()
0.58736140513792634
>>> ud2 = galsim.UniformDeviate(215324)
>>> ud2()
0.58736140513792634
"""
def __init__(self, seed=None):
self._rng_type = _galsim.BaseDeviateImpl
self._rng_args = ()
self.reset(seed)
[docs] def seed(self, seed=0):
"""Seed the pseudo-random number generator with a given integer value.
Parameters:
seed: An int value to be used to seed the random number generator. Using 0
means to generate a seed from the system. [default: 0]
"""
if seed == int(seed):
self._seed(int(seed))
else:
raise TypeError("BaseDeviate seed must be an integer. Got %s"%seed)
[docs] def _seed(self, seed=0):
"""Equivalent to `seed`, but without any type checking.
"""
self._rng.seed(seed)
[docs] def reset(self, seed=None):
"""Reset the pseudo-random number generator, severing connections to any other deviates.
Providing another `BaseDeviate` object as the seed connects this deviate with the other
one, so they will both use the same underlying random number generator.
Parameters:
seed: Something that can seed a `BaseDeviate`: an integer seed or another
`BaseDeviate`. Using None means to generate a seed from the system.
[default: None]
"""
if isinstance(seed, BaseDeviate):
self._reset(seed)
elif isinstance(seed, str):
self._rng = self._rng_type(_galsim.BaseDeviateImpl(seed), *self._rng_args)
elif seed is None:
self._rng = self._rng_type(_galsim.BaseDeviateImpl(0), *self._rng_args)
elif isinteger(seed):
self._rng = self._rng_type(_galsim.BaseDeviateImpl(int(seed)), *self._rng_args)
else:
raise TypeError("BaseDeviate must be initialized with either an int or another "
"BaseDeviate")
[docs] def _reset(self, rng):
"""Equivalent to `reset`, but rng must be a `BaseDeviate` (not an int), and there
is no type checking.
"""
self._rng = self._rng_type(rng._rng, *self._rng_args)
@property
def np(self):
"""Shorthand for self.as_numpy_generator()
"""
return self.as_numpy_generator()
[docs] def as_numpy_generator(self):
"""Return a numpy.random.Generator object that uses the current BaseDeviate for the
underlying bit generations.
This allows you to use the (probably) more familiar numpy functions, while maintaining
GalSim's guarantees about random number stability across platforms.
Example::
>>> rng = galsim.BaseDeviate(1234)
>>> gen = rng.as_numpy_generator()
>>> uniform = gen.uniform(1, 10, size=10)
>>> norm = gen.normal(0, 3, size=20)
There is also a shorthand syntax that may be convenient.
The property `np` is equivalent to this method, so you can also write::
>>> uniform = rng.np.uniform(1, 10, size=10)
>>> norm = rng.np.normal(0, 3, size=20)
"""
return np.random.Generator(GalSimBitGenerator(self))
[docs] def duplicate(self):
"""Create a duplicate of the current `BaseDeviate` object.
The subsequent series from each copy of the `BaseDeviate` will produce identical values::
>>> u = galsim.UniformDeviate(31415926)
>>> u()
0.17100770119577646
>>> u2 = u.duplicate()
>>> u()
0.49095047544687986
>>> u()
0.10306670609861612
>>> u2()
0.49095047544687986
>>> u2()
0.10306670609861612
>>> u2()
0.13129289541393518
>>> u()
0.13129289541393518
"""
ret = BaseDeviate.__new__(self.__class__)
ret.__dict__.update(self.__dict__)
ret._rng = self._rng.duplicate()
return ret
def __copy__(self):
return self.duplicate()
def __getstate__(self):
d = self.__dict__.copy()
d['rng_str'] = self.serialize()
d.pop('_rng')
return d
def __setstate__(self, d):
self.__dict__ = d
rng = _galsim.BaseDeviateImpl(d['rng_str'])
self._rng = self._rng_type(rng, *self._rng_args)
[docs] def clearCache(self):
"""Clear the internal cache of the `BaseDeviate`, if any. This is currently only relevant
for `GaussianDeviate`, since it generates two values at a time, saving the second one to
use for the next output value.
"""
self._rng.clearCache()
[docs] def discard(self, n, suppress_warnings=False):
"""Discard n values from the current sequence of pseudo-random numbers.
This is typically used to keep two random number generators in sync when one of them
is used to generate some random values. The other can quickly discard the same number
of random values to stay in sync with the first one.
Parameters:
n: The number of raw random numbers to discard.
suppress_warnings: Whether to suppress warnings related to detecting when this
action is not likely to reliably keep two random number
generators in sync. [default: False]
"""
if not self.has_reliable_discard and not suppress_warnings:
galsim_warn(self.__class__.__name__ +
" does not use a consistent number of randoms per generated value, " +
"so discard cannot be guaranteed to keep two random deviates in sync.")
if n%2 == 1 and self.generates_in_pairs and not suppress_warnings:
galsim_warn(self.__class__.__name__ +
" uses two randoms per pair of generated values, so discarding " +
"an odd number of randoms probably doesn't make sense.")
self._rng.discard(int(n))
@property
def has_reliable_discard(self):
"""Indicates whether the generator always uses 1 rng per value.
If it does, then `discard` can reliably be used to keep two generators in sync when one
of them generated some values and the other didn't.
This is False for PoissonDeviate, Chi2Deviate, and GammaDeviate.
See also `BaseDeviate.generates_in_pairs`.
"""
return True
@property
def generates_in_pairs(self):
"""Indicates whether the generator uses 2 rngs values per 2 returned values.
GaussianDeviate has a slight wrinkle to the `BaseDeviate.has_reliable_discard` story.
It always uses two rng values to generate two Gaussian deviates. So if the number of
generated values is even, then discard can keep things in sync. However, if an odd
number of values are generated, then you to generate one more value (and discard it)
to keep things synced up.
This is only True for GaussianDeviate.
"""
return False
[docs] def raw(self):
"""Generate the next pseudo-random number and rather than return the appropriate kind
of random deviate for this class, just return the raw integer value that would have been
used to generate this value.
"""
return int(self._rng.raw())
[docs] def generate(self, array):
"""Generate many pseudo-random values, filling in the values of a numpy array.
"""
array_1d = np.ascontiguousarray(array.ravel(),dtype=float)
#assert(array_1d.strides[0] == array_1d.itemsize)
_a = array_1d.__array_interface__['data'][0]
self._rng.generate(len(array_1d), _a)
if array_1d.data != array.data:
# array_1d is not a view into the original array. Need to copy back.
np.copyto(array, array_1d.reshape(array.shape), casting='unsafe')
[docs] def add_generate(self, array):
"""Generate many pseudo-random values, adding them to the values of a numpy array.
"""
array_1d = np.ascontiguousarray(array.ravel(),dtype=float)
#assert(array_1d.strides[0] == array_1d.itemsize)
_a = array_1d.__array_interface__['data'][0]
self._rng.add_generate(len(array_1d), _a)
if array_1d.data != array.data:
# array_1d is not a view into the original array. Need to copy back.
np.copyto(array, array_1d.reshape(array.shape), casting='unsafe')
def __eq__(self, other):
return (self is other or
(isinstance(other, self.__class__) and
self._rng_type == other._rng_type and
self._rng_args == other._rng_args and
self.serialize() == other.serialize()))
def __ne__(self, other):
return not self.__eq__(other)
__hash__ = None
def serialize(self):
return str(self._rng.serialize())
def _seed_repr(self):
s = self.serialize().split(' ')
return " ".join(s[:3])+" ... "+" ".join(s[-3:])
def __repr__(self):
return "galsim.BaseDeviate(%r)"%self._seed_repr()
def __str__(self):
return "galsim.BaseDeviate(%r)"%self._seed_repr()
[docs]class GaussianDeviate(BaseDeviate):
"""Pseudo-random number generator with Gaussian distribution.
See http://en.wikipedia.org/wiki/Gaussian_distribution for further details.
Successive calls to ``g()`` generate pseudo-random values distributed according to a Gaussian
distribution with the provided ``mean``, ``sigma``::
>>> g = galsim.GaussianDeviate(31415926)
>>> g()
0.5533754000847082
>>> g()
1.0218588970190354
Parameters:
seed: Something that can seed a `BaseDeviate`: an integer seed or another
`BaseDeviate`. Using 0 means to generate a seed from the system.
[default: None]
mean: Mean of Gaussian distribution. [default: 0.]
sigma: Sigma of Gaussian distribution. [default: 1.; Must be > 0]
"""
def __init__(self, seed=None, mean=0., sigma=1.):
if sigma < 0.:
raise GalSimRangeError("GaussianDeviate sigma must be > 0.", sigma, 0.)
self._rng_type = _galsim.GaussianDeviateImpl
self._rng_args = (float(mean), float(sigma))
self.reset(seed)
@property
def mean(self):
"""The mean of the Gaussian distribution.
"""
return self._rng_args[0]
@property
def sigma(self):
"""The sigma of the Gaussian distribution.
"""
return self._rng_args[1]
@property
def generates_in_pairs(self):
return True
[docs] def __call__(self):
"""Draw a new random number from the distribution.
Returns a Gaussian deviate with the given mean and sigma.
"""
return self._rng.generate1()
[docs] def generate_from_variance(self, array):
"""Generate many Gaussian deviate values using the existing array values as the
variance for each.
"""
array_1d = np.ascontiguousarray(array.ravel(), dtype=float)
#assert(array_1d.strides[0] == array_1d.itemsize)
_a = array_1d.__array_interface__['data'][0]
self._rng.generate_from_variance(len(array_1d), _a)
if array_1d.data != array.data:
# array_1d is not a view into the original array. Need to copy back.
np.copyto(array, array_1d.reshape(array.shape), casting='unsafe')
def __repr__(self):
return 'galsim.GaussianDeviate(seed=%r, mean=%r, sigma=%r)'%(
self._seed_repr(), self.mean, self.sigma)
def __str__(self):
return 'galsim.GaussianDeviate(mean=%r, sigma=%r)'%(self.mean, self.sigma)
[docs]class BinomialDeviate(BaseDeviate):
"""Pseudo-random Binomial deviate for ``N`` trials each of probability ``p``.
``N`` is number of 'coin flips,' ``p`` is probability of 'heads,' and each call returns an
integer value where 0 <= value <= N gives the number of heads. See
http://en.wikipedia.org/wiki/Binomial_distribution for more information.
Successive calls to ``b()`` generate pseudo-random integer values distributed according to a
binomial distribution with the provided ``N``, ``p``::
>>> b = galsim.BinomialDeviate(31415926, N=10, p=0.3)
>>> b()
2
>>> b()
3
Parameters:
seed: Something that can seed a `BaseDeviate`: an integer seed or another
`BaseDeviate`. Using 0 means to generate a seed from the system.
[default: None]
N: The number of 'coin flips' per trial. [default: 1; Must be > 0]
p: The probability of success per coin flip. [default: 0.5; Must be > 0]
"""
def __init__(self, seed=None, N=1, p=0.5):
self._rng_type = _galsim.BinomialDeviateImpl
self._rng_args = (int(N), float(p))
self.reset(seed)
@property
def n(self):
"""The number of 'coin flips'.
"""
return self._rng_args[0]
@property
def p(self):
"""The probability of success per 'coin flip'.
"""
return self._rng_args[1]
[docs] def __call__(self):
"""Draw a new random number from the distribution.
Returns a Binomial deviate with the given n and p.
"""
return self._rng.generate1()
def __repr__(self):
return 'galsim.BinomialDeviate(seed=%r, N=%r, p=%r)'%(self._seed_repr(), self.n, self.p)
def __str__(self):
return 'galsim.BinomialDeviate(N=%r, p=%r)'%(self.n, self.p)
[docs]class PoissonDeviate(BaseDeviate):
"""Pseudo-random Poisson deviate with specified ``mean``.
The input ``mean`` sets the mean and variance of the Poisson deviate. An integer deviate with
this distribution is returned after each call.
See http://en.wikipedia.org/wiki/Poisson_distribution for more details.
Successive calls to ``p()`` generate pseudo-random integer values distributed according to a
Poisson distribution with the specified ``mean``::
>>> p = galsim.PoissonDeviate(31415926, mean=100)
>>> p()
94
>>> p()
106
Parameters:
seed: Something that can seed a `BaseDeviate`: an integer seed or another
`BaseDeviate`. Using 0 means to generate a seed from the system.
[default: None]
mean: Mean of the distribution. [default: 1; Must be > 0]
"""
def __init__(self, seed=None, mean=1.):
if mean < 0:
raise GalSimValueError("PoissonDeviate is only defined for mean >= 0.", mean)
self._rng_type = _galsim.PoissonDeviateImpl
self._rng_args = (float(mean),)
self.reset(seed)
@property
def mean(self):
"""The mean of the distribution.
"""
return self._rng_args[0]
@property
def has_reliable_discard(self):
return False
[docs] def __call__(self):
"""Draw a new random number from the distribution.
Returns a Poisson deviate with the given mean.
"""
return self._rng.generate1()
[docs] def generate_from_expectation(self, array):
"""Generate many Poisson deviate values using the existing array values as the
expectation value (aka mean) for each.
"""
if np.any(array < 0):
raise GalSimValueError("Expectation array may not have values < 0.", array)
array_1d = np.ascontiguousarray(array.ravel(), dtype=float)
#assert(array_1d.strides[0] == array_1d.itemsize)
_a = array_1d.__array_interface__['data'][0]
self._rng.generate_from_expectation(len(array_1d), _a)
if array_1d.data != array.data:
# array_1d is not a view into the original array. Need to copy back.
np.copyto(array, array_1d.reshape(array.shape), casting='unsafe')
def __repr__(self):
return 'galsim.PoissonDeviate(seed=%r, mean=%r)'%(self._seed_repr(), self.mean)
def __str__(self):
return 'galsim.PoissonDeviate(mean=%r)'%(self.mean)
[docs]class WeibullDeviate(BaseDeviate):
"""Pseudo-random Weibull-distributed deviate for shape parameter ``a`` and scale parameter ``b``.
The Weibull distribution is related to a number of other probability distributions; in
particular, it interpolates between the exponential distribution (a=1) and the Rayleigh
distribution (a=2).
See http://en.wikipedia.org/wiki/Weibull_distribution (a=k and b=lambda in the notation adopted
in the Wikipedia article) for more details. The Weibull distribution is real valued and
produces deviates >= 0.
Successive calls to ``w()`` generate pseudo-random values distributed according to a Weibull
distribution with the specified shape and scale parameters ``a`` and ``b``::
>>> w = galsim.WeibullDeviate(31415926, a=1.3, b=4)
>>> w()
1.1038481241018219
>>> w()
2.957052966368049
Parameters:
seed: Something that can seed a `BaseDeviate`: an integer seed or another
`BaseDeviate`. Using 0 means to generate a seed from the system.
[default: None]
a: Shape parameter of the distribution. [default: 1; Must be > 0]
b: Scale parameter of the distribution. [default: 1; Must be > 0]
"""
def __init__(self, seed=None, a=1., b=1.):
self._rng_type = _galsim.WeibullDeviateImpl
self._rng_args = (float(a), float(b))
self.reset(seed)
@property
def a(self):
"""The shape parameter, a.
"""
return self._rng_args[0]
@property
def b(self):
"""The scale parameter, b.
"""
return self._rng_args[1]
[docs] def __call__(self):
"""Draw a new random number from the distribution.
Returns a Weibull-distributed deviate with the given shape parameters a and b.
"""
return self._rng.generate1()
def __repr__(self):
return 'galsim.WeibullDeviate(seed=%r, a=%r, b=%r)'%(self._seed_repr(), self.a, self.b)
def __str__(self):
return 'galsim.WeibullDeviate(a=%r, b=%r)'%(self.a, self.b)
[docs]class GammaDeviate(BaseDeviate):
"""A Gamma-distributed deviate with shape parameter ``k`` and scale parameter ``theta``.
See http://en.wikipedia.org/wiki/Gamma_distribution.
(Note: we use the k, theta notation. If you prefer alpha, beta, use k=alpha, theta=1/beta.)
The Gamma distribution is a real valued distribution producing deviates >= 0.
Successive calls to ``g()`` generate pseudo-random values distributed according to a gamma
distribution with the specified shape and scale parameters ``k`` and ``theta``::
>>> gam = galsim.GammaDeviate(31415926, k=1, theta=2)
>>> gam()
0.37508882726316
>>> gam()
1.3504199388358704
Parameters:
seed: Something that can seed a `BaseDeviate`: an integer seed or another
`BaseDeviate`. Using 0 means to generate a seed from the system.
[default: None]
k: Shape parameter of the distribution. [default: 1; Must be > 0]
theta: Scale parameter of the distribution. [default: 1; Must be > 0]
"""
def __init__(self, seed=None, k=1., theta=1.):
self._rng_type = _galsim.GammaDeviateImpl
self._rng_args = (float(k), float(theta))
self.reset(seed)
@property
def k(self):
"""The shape parameter, k.
"""
return self._rng_args[0]
@property
def theta(self):
"""The scale parameter, theta.
"""
return self._rng_args[1]
@property
def has_reliable_discard(self):
return False
[docs] def __call__(self):
"""Draw a new random number from the distribution.
Returns a Gamma-distributed deviate with the given k and theta.
"""
return self._rng.generate1()
def __repr__(self):
return 'galsim.GammaDeviate(seed=%r, k=%r, theta=%r)'%(
self._seed_repr(), self.k, self.theta)
def __str__(self):
return 'galsim.GammaDeviate(k=%r, theta=%r)'%(self.k, self.theta)
[docs]class Chi2Deviate(BaseDeviate):
"""Pseudo-random Chi^2-distributed deviate for degrees-of-freedom parameter ``n``.
See http://en.wikipedia.org/wiki/Chi-squared_distribution (note that k=n in the notation
adopted in the Boost.Random routine called by this class). The Chi^2 distribution is a
real-valued distribution producing deviates >= 0.
Successive calls to ``chi2()`` generate pseudo-random values distributed according to a
chi-square distribution with the specified degrees of freedom, ``n``::
>>> chi2 = galsim.Chi2Deviate(31415926, n=7)
>>> chi2()
7.9182211987712385
>>> chi2()
6.644121724269535
Parameters:
seed: Something that can seed a `BaseDeviate`: an integer seed or another
`BaseDeviate`. Using 0 means to generate a seed from the system.
[default: None]
n: Number of degrees of freedom for the output distribution. [default: 1;
Must be > 0]
"""
def __init__(self, seed=None, n=1.):
self._rng_type = _galsim.Chi2DeviateImpl
self._rng_args = (float(n),)
self.reset(seed)
@property
def n(self):
"""The number of degrees of freedom.
"""
return self._rng_args[0]
@property
def has_reliable_discard(self):
return False
[docs] def __call__(self):
"""Draw a new random number from the distribution.
Returns a Chi2-distributed deviate with the given number of degrees of freedom.
"""
return self._rng.generate1()
def __repr__(self):
return 'galsim.Chi2Deviate(seed=%r, n=%r)'%(self._seed_repr(), self.n)
def __str__(self):
return 'galsim.Chi2Deviate(n=%r)'%(self.n)
[docs]class DistDeviate(BaseDeviate):
"""A class to draw random numbers from a user-defined probability distribution.
DistDeviate is a `BaseDeviate` class that can be used to draw from an arbitrary probability
distribution. The probability distribution passed to DistDeviate can be given one of three
ways: as the name of a file containing a 2d ASCII array of x and P(x), as a `LookupTable`
mapping x to P(x), or as a callable function.
Once given a probability, DistDeviate creates a table of the cumulative probability and draws
from it using a `UniformDeviate`. The precision of its outputs can be controlled with the
keyword ``npoints``, which sets the number of points DistDeviate creates for its internal table
of CDF(x). To prevent errors due to non-monotonicity, the interpolant for this internal table
is always linear.
Two keywords, ``x_min`` and ``x_max``, define the support of the function. They must be passed
if a callable function is given to DistDeviate, unless the function is a `LookupTable`, which
has its own defined endpoints. If a filename or `LookupTable` is passed to DistDeviate, the
use of ``x_min`` or ``x_max`` will result in an error.
If given a table in a file, DistDeviate will construct an interpolated `LookupTable` to obtain
more finely gridded probabilities for generating the cumulative probability table. The default
``interpolant`` is linear, but any interpolant understood by `LookupTable` may be used. We
caution against the use of splines because they can cause non-monotonic behavior. Passing the
``interpolant`` keyword next to anything but a table in a file will result in an error.
**Examples**:
Some sample initialization calls::
>>> d = galsim.DistDeviate(function=f, x_min=x_min, x_max=x_max)
Initializes d to be a DistDeviate instance with a distribution given by the callable function
``f(x)`` from ``x=x_min`` to ``x=x_max`` and seeds the PRNG using current time::
>>> d = galsim.DistDeviate(1062533, function=file_name, interpolant='floor')
Initializes d to be a DistDeviate instance with a distribution given by the data in file
``file_name``, which must be a 2-column ASCII table, and seeds the PRNG using the integer
seed 1062533. It generates probabilities from ``file_name`` using the interpolant 'floor'::
>>> d = galsim.DistDeviate(rng, function=galsim.LookupTable(x,p))
Initializes d to be a DistDeviate instance with a distribution given by P(x), defined as two
arrays ``x`` and ``p`` which are used to make a callable `LookupTable`, and links the
DistDeviate PRNG to the already-existing random number generator ``rng``.
Successive calls to ``d()`` generate pseudo-random values with the given probability
distribution::
>>> d = galsim.DistDeviate(31415926, function=lambda x: 1-abs(x), x_min=-1, x_max=1)
>>> d()
-0.4151921102709466
>>> d()
-0.00909781188974034
Parameters:
seed: Something that can seed a `BaseDeviate`: an integer seed or another
`BaseDeviate`. Using 0 means to generate a seed from the system.
[default: None]
function: A callable function giving a probability distribution or the name of a
file containing a probability distribution as a 2-column ASCII table.
[required]
x_min: The minimum desired return value (required for non-`LookupTable`
callable functions; will raise an error if not passed in that case, or if
passed in any other case) [default: None]
x_max: The maximum desired return value (required for non-`LookupTable`
callable functions; will raise an error if not passed in that case, or if
passed in any other case) [default: None]
interpolant: Type of interpolation used for interpolating a file (causes an error if
passed alongside a callable function). Options are given in the
documentation for `LookupTable`. [default: 'linear']
npoints: Number of points DistDeviate should create for its internal interpolation
tables. [default: 256, unless the function is a non-log `LookupTable`, in
which case it uses the table's x values]
clip_neg: Clip any negative input values to zero. [default: False; an error will
be raised if any negative probabilities are found.]
"""
def __init__(self, seed=None, function=None, x_min=None,
x_max=None, interpolant=None, npoints=None, clip_neg=False):
# Set up the PRNG
self._rng_type = _galsim.UniformDeviateImpl
self._rng_args = ()
self.reset(seed)
# Basic input checking and setups
if function is None:
raise TypeError('You must pass a function to DistDeviate!')
self._interpolant = interpolant
self._npoints = npoints
self._xmin = x_min
self._xmax = x_max
# Figure out if a string is a filename or something we should be using in an eval call
if isinstance(function, str):
self._function = function # Save the inputs to be used in repr
if os.path.isfile(function):
if interpolant is None:
interpolant='linear'
if x_min or x_max:
raise GalSimIncompatibleValuesError(
"Cannot pass x_min or x_max with a filename argument",
function=function, x_min=x_min, x_max=x_max)
function = LookupTable.from_file(function, interpolant=interpolant)
x_min = function.x_min
x_max = function.x_max
else:
try:
function = math_eval('lambda x : ' + function)
if x_min is not None: # is not None in case x_min=0.
function(x_min)
else:
# Somebody would be silly to pass a string for evaluation without x_min,
# but we'd like to throw reasonable errors in that case anyway
function(0.6) # A value unlikely to be a singular point of a function
except Exception as e:
raise GalSimValueError(
"String function must either be a valid filename or something that "
"can eval to a function of x.\n"
"Caught error: {0}".format(e), self._function)
else:
# Check that the function is actually a function
if not hasattr(function, '__call__'):
raise TypeError('function must be a callable function or a string')
if interpolant:
raise GalSimIncompatibleValuesError(
"Cannot provide an interpolant with a callable function argument",
interpolant=interpolant, function=function)
if isinstance(function, LookupTable):
if (x_min not in (None, function.x_min)) or (x_max not in (None, function.x_max)):
raise GalSimIncompatibleValuesError(
"Cannot provide x_min or x_max with a LookupTable function",
function=function, x_min=x_min, x_max=x_max)
x_min = function.x_min
x_max = function.x_max
else:
if x_min is None or x_max is None:
raise GalSimIncompatibleValuesError(
"Must provide x_min and x_max when function argument is a regular "
"python callable function",
function=function, x_min=x_min, x_max=x_max)
self._function = function # Save the inputs to be used in repr
# Compute the probability distribution function, pdf(x)
if (npoints is None and isinstance(function, LookupTable) and
not function.x_log and not function.f_log):
xarray = np.array(function.x, dtype=float)
pdf = np.array(function.f, dtype=float)
# Set up pdf, so cumsum basically does a cumulative trapz integral
# On Python 3.4, doing pdf[1:] += pdf[:-1] the last value gets messed up.
# Writing it this way works. (Maybe slightly slower though, so if we stop
# supporting python 3.4, consider switching to the += version.)
pdf[1:] = pdf[1:] + pdf[:-1]
pdf[1:] *= np.diff(xarray)
pdf[0] = 0.
else:
if npoints is None: npoints = 256
xarray = x_min+(1.*x_max-x_min)/(npoints-1)*np.array(range(npoints),float)
# Integrate over the range of x in case the function is doing something weird here.
pdf = [0.] + [integ.int1d(function, xarray[i], xarray[i+1])
for i in range(npoints - 1)]
pdf = np.array(pdf)
# Check that the probability is nonnegative
if clip_neg:
# Write it this way so nan -> 0 as well as negative values.
w = np.where(~(pdf >= 0))
pdf[w] = 0.
elif not np.all(pdf >= 0.):
raise GalSimValueError('Negative probability found in DistDeviate.',function)
# Compute the cumulative distribution function = int(pdf(x),x)
cdf = np.cumsum(pdf)
# Quietly renormalize the probability if it wasn't already normalized
totalprobability = cdf[-1]
cdf /= totalprobability
self._inverse_cdf = LookupTable(cdf, xarray, interpolant='linear')
self.x_min = x_min
self.x_max = x_max
[docs] def val(self, p):
r"""
Return the value :math:`x` of the input function to `DistDeviate` such that ``p`` =
:math:`F(x)`, where :math:`F` is the cumulattive probability distribution function:
.. math::
F(x) = \int_{-\infty}^x \mathrm{pdf}(t) dt
This function is typically called by `__call__`, which generates a random p
between 0 and 1 and calls ``self.val(p)``.
Parameters:
p: The desired cumulative probabilty p.
Returns:
the corresponding x such that :math:`p = F(x)`.
"""
if p<0 or p>1:
raise GalSimRangeError('Invalid cumulative probability for DistDeviate', p, 0., 1.)
return self._inverse_cdf(p)
[docs] def __call__(self):
"""Draw a new random number from the distribution.
"""
return self._inverse_cdf(self._rng.generate1())
[docs] def generate(self, array):
"""Generate many pseudo-random values, filling in the values of a numpy array.
"""
p = np.empty_like(array)
BaseDeviate.generate(self, p) # Fill with unform deviate values
np.copyto(array, self._inverse_cdf(p)) # Convert from p -> x
[docs] def add_generate(self, array):
"""Generate many pseudo-random values, adding them to the values of a numpy array.
"""
p = np.empty_like(array)
BaseDeviate.generate(self, p)
array += self._inverse_cdf(p)
def __repr__(self):
return ('galsim.DistDeviate(seed=%r, function=%r, x_min=%r, x_max=%r, interpolant=%r, '
'npoints=%r)')%(self._seed_repr(), self._function, self._xmin, self._xmax,
self._interpolant, self._npoints)
def __str__(self):
return 'galsim.DistDeviate(function="%s", x_min=%s, x_max=%s, interpolant=%s, npoints=%s)'%(
self._function, self._xmin, self._xmax, self._interpolant, self._npoints)
def __eq__(self, other):
return (self is other or
(isinstance(other, DistDeviate) and
self.serialize() == other.serialize() and
self._function == other._function and
self._xmin == other._xmin and
self._xmax == other._xmax and
self._interpolant == other._interpolant and
self._npoints == other._npoints))
class GalSimBitGenerator(np.random.BitGenerator):
"""A numpy.random.BitGenerator that uses the GalSim C++-layer random number generator
for the random bit generation.
Parameters:
rng: The galsim.BaseDeviate object to use for the underlying bit generation.
"""
def __init__(self, rng):
super().__init__(0)
self.rng = rng
self.rng._rng.setup_bitgen(self.capsule)
def permute(rng, *args):
"""Randomly permute one or more lists.
If more than one list is given, then all lists will have the same random permutation
applied to it.
Parameters:
rng: The random number generator to use. (This will be converted to a `UniformDeviate`.)
args: Any number of lists to be permuted.
"""
ud = UniformDeviate(rng)
if len(args) == 0:
raise TypeError("permute called with no lists to permute")
# We use an algorithm called the Knuth shuffle, which is based on the Fisher-Yates shuffle.
# See http://en.wikipedia.org/wiki/Fisher-Yates_shuffle for more information.
n = len(args[0])
for i in range(n-1,1,-1):
j = int((i+1) * ud())
if j == i+1: j = i # I'm not sure if this is possible, but just in case...
for lst in args:
lst[i], lst[j] = lst[j], lst[i]