Lookup Tables

class galsim.LookupTable(x, f, interpolant='spline', x_log=False, f_log=False)[source]

LookupTable represents a lookup table to store function values that may be slow to calculate, for which interpolating from a lookup table is sufficiently accurate.

A LookupTable may be constructed from two arrays (lists, tuples, or NumPy arrays of floats/doubles):

>>> args = [...]
>>> vals = []
>>> for arg in args:
...     val = calculateVal(arg)
...     vals.append(val)
>>> table = galsim.LookupTable(x=args,f=vals)

Then you can use this table as a replacement for the slow calculation:

>>> other_args = [...]
>>> for arg in other_args:
...     val = table(arg)
...     [... use val ...]

The default interpolation method is a natural cubic spline. This is usually the best choice, but we also provide other options, which can be specified by the interpolant kwarg. The choices include ‘floor’, ‘ceil’, ‘linear’, ‘spline’, or a galsim.Interpolant object:

  • ‘floor’ takes the value from the previous argument in the table.

  • ‘ceil’ takes the value from the next argument in the table.

  • ‘nearest’ takes the value from the nearest argument in the table.

  • ‘linear’ does linear interpolation between these two values.

  • ‘spline’ uses a cubic spline interpolation, so the interpolated values are smooth at each argument in the table.

  • a galsim.Interpolant object or a string convertible to one. This option can be used for Lanczos or Quintic interpolation, for example.

Note that specifying the string ‘nearest’ or ‘linear’ will use a LookupTable-optimized interpolant instead of galsim.Nearest or galsim.Linear, though the latter options can still be used by passing an Interpolant object instead of a string. Also note that to use a galsim.Interpolant in a LookupTable, the input data must be equally spaced, or logarithmically spaced if x_log is set to True (see below). Finally, although natural cubic spline used when interpolant=’spline’ and the cubic convolution interpolant used when the interpolant is galsim.Cubic both produce piecewise cubic polynomial interpolations, their treatments of the continuity of derivatives are different (the natural spline is smoother).

There are also two factory functions which can be used to build a LookupTable:

LookupTable.from_func

makes a LookupTable from a callable function

LookupTable.from_file

reads in a file of x and f values.

The user can also opt to interpolate in log(x) and/or log(f) (if not using a galsim.Interpolant), though this is not the default. It may be a wise choice depending on the particular function, e.g., for a nearly power-law f(x) (or at least one that is locally power-law-ish for much of the x range) then it might be a good idea to interpolate in log(x) and log(f) rather than x and f.

Parameters:

  • x – The list, tuple, or NumPy array of x values.

  • f – The list, tuple, or NumPy array of f(x) values.

  • interpolant – Type of interpolation to use, with the options being ‘floor’, ‘ceil’, ‘nearest’, ‘linear’, ‘spline’, or a galsim.Interpolant or string convertible to one. [default: ‘spline’]

  • x_log – Set to True if you wish to interpolate using log(x) rather than x. Note that all inputs / outputs will still be x, it’s just a question of how the interpolation is done. [default: False]

  • f_log – Set to True if you wish to interpolate using log(f) rather than f. Note that all inputs / outputs will still be f, it’s just a question of how the interpolation is done. [default: False]

__call__(x)[source]

Interpolate the LookupTable to get f(x) at some x value(s).

When the LookupTable object is called with a single argument, it returns the value at that argument. An exception will be thrown automatically if the x value is outside the range of the original tabulated values. The value that is returned is the same type as that provided as an argument, e.g., if a single value x is provided then a single value of f is returned; if a tuple of x values is provided then a tuple of f values is returned; and so on. Even if interpolation was done using the x_log option, the user should still provide x rather than log(x).

Parameters:

x – The x value(s) for which f(x) should be calculated via interpolation on the original (x,f) lookup table. x can be a single float/double, or a tuple, list, or arbitrarily shaped 1- or 2-dimensional NumPy array.

Returns:

the interpolated f(x) value(s).

classmethod from_file(file_name, interpolant='spline', x_log=False, f_log=False, amplitude=1.0)[source]

Create a LookupTable from a file of x, f values.

This reads in a file, which should contain two columns with the x and f values.

Parameters:

  • file_name – A file from which to read the (x,f) pairs.

  • interpolant – Type of interpolation to use. [default: ‘spline’]

  • x_log – Whether the x values should be uniform in log rather than lienar. [default: False]

  • f_log – Whether the f values should be interpolated using their logarithms rather than their raw values. [default: False]

  • amplitude – An optional scaling of the f values relative to the values in the file [default: 1.0]

classmethod from_func(func, x_min, x_max, npoints=2000, interpolant='spline', x_log=False, f_log=False)[source]

Create a LookupTable from a callable function

This constructs a LookupTable over the given range from x_min and x_max, calculating the corresponding f values from the given function (technically any callable object).

Parameters:

  • func – A callable function.

  • x_min – The minimum x value at which to evalue the function and store in the lookup table.

  • x_max – The maximum x value at which to evalue the function and store in the lookup table.

  • npoints – Number of x values at which to evaluate the function. [default: 2000]

  • interpolant – Type of interpolation to use. [default: ‘spline’]

  • x_log – Whether the x values should be uniform in log rather than lienar. [default: False]

  • f_log – Whether the f values should be interpolated using their logarithms rather than their raw values. [default: False]

integrate(x_min=None, x_max=None)[source]

Calculate an estimate of the integral of the tabulated function from x_min to x_max:

\[\int_{x_\mathrm{min}}^{x_\mathrm{max}} f(x) dx\]

This function is not implemented for LookupTables that use log for either x or f, or that use a galsim.Interpolant. Also, if x_min or x_max are beyond the range of the tabulated function, the function will be considered to be zero there.

Note

The simplest version of this function is equivalent in functionality to the numpy trapz function. However, it is usually significantly faster. If you have a time-critical integration for which you are currently using np.trapz:

>>> ans = np.trapz(f, x)

the following replacement may be faster:

>>> ans = galsim.trapz(f, x)

which is an alias for:

>>> ans = galsim._LookupTable(x, f, 'linear').integrate()
Parameters:

  • x_min – The minimum abscissa to use for the integral. [default: None, which means to use self.x_min]

  • x_max – The maximum abscissa to use for the integral. [default: None, which means to use self.x_max]

Returns:

an estimate of the integral

integrate_product(g, x_min=None, x_max=None, x_factor=1.0)[source]

Calculate an estimate of the integral of the tabulated function multiplied by a second function from x_min to x_max:

\[\int_{x_\mathrm{min}}^{x_\mathrm{max}} f(x) g(x) dx\]

If the second function, \(g(x)\), is another LookupTable, then the quadrature will use the abscissae from both that function and \(f(x)\) (i.e. self). Otherwise, the second function will be evaluated at the abscissae of \(f(x)\).

This function is not implemented for LookupTables that use log for either x or f, or that use a galsim.Interpolant. Also, if x_min or x_max are beyond the range of either tabulated function, the function will be considered to be zero there.

Also, the second function \(g(x)\) is always approximated with linear interpolation between the abscissae, even if it is a LookupTable with a different specified interpolation.

Parameters:

  • g – The function to multiply by the current function for the integral.

  • x_min – The minimum abscissa to use for the integral. [default: None, which means to use self.x_min]

  • x_max – The maximum abscissa to use for the integral. [default: None, which means to use self.x_max]

  • x_factor – Optionally scale the x values of f by this factor when doing the integral. I.e. Find \(\int f(x x_\mathrm{factor}) g(x) dx\). [default: 1]

Returns:

an estimate of the integral

property x_max

The maximum x value in the lookup table.

property x_min

The minimum x value in the lookup table.

class galsim.LookupTable2D(x, y, f, dfdx=None, dfdy=None, d2fdxdy=None, interpolant='linear', edge_mode='raise', constant=0)[source]

LookupTable2D represents a 2-dimensional lookup table to store function values that may be slow to calculate, for which interpolating from a lookup table is sufficiently accurate. A LookupTable2D is also useful for evaluating periodic 2-d functions given samples from a single period.

A LookupTable2D representing the function f(x, y) may be constructed from a list or array of x values, a list or array of y values, and a 2D array of function evaluations at all combinations of x and y values. For instance:

>>> x = np.arange(5)
>>> y = np.arange(8)
>>> z = x + y[:, np.newaxis]  # function is x + y, dimensions of z are (8, 5)
>>> tab2d = galsim.LookupTable2D(x, y, z)

To evaluate new function values with the lookup table, use the () operator:

>>> print tab2d(2.2, 3.3)
5.5

The () operator can also accept sequences (lists, tuples, numpy arrays, …) for the x and y arguments at which to evaluate the LookupTable2D. Normally, the x and y sequences should have the same length, which will also be the length of the output sequence:

>>> print tab2d([1, 2], [3, 5])
array([ 4., 7.])

If you add grid=True as an additional kwarg, however, then the () operator will generate interpolated values at the outer product of x-values and y-values. So in this case, the x and y sequences can have different lengths Nx and Ny, and the result will be a 2D array with dimensions (Nx, Ny):

>>> print tab2d([1, 2], [3, 5], grid=True)
array([[ 4., 6.],
       [ 5., 7.]])

The default interpolation method is linear. Other choices for the interpolant are:

  • ‘floor’

  • ‘ceil’

  • ‘nearest’

  • ‘spline’ (a Catmull-Rom cubic spline).

  • a galsim.Interpolant or string convertible to one.

>>> tab2d = galsim.LookupTable2D(x, y, z, interpolant='floor')
>>> tab2d(2.2, 3.7)
5.0
>>> tab2d = galsim.LookupTable2D(x, y, z, interpolant='ceil')
>>> tab2d(2.2, 3.7)
7.0
>>> tab2d = galsim.LookupTable2D(x, y, z, interpolant='nearest')
>>> tab2d(2.2, 3.7)
6.0

For interpolant=’spline’ or a galsim.Interpolant, the input arrays must be uniformly spaced. For interpolant=’spline’, the derivatives df / dx, df / dy, and d^2 f / dx dy at grid-points may also optionally be provided if they’re known, which will generally yield a more accurate interpolation (these derivatives will be estimated from finite differences if they’re not provided).

The edge_mode keyword describes how to handle extrapolation beyond the initial input range. Possibilities include:

  • ‘raise’: raise an exception. (This is the default.)

  • ‘warn’: issues a warning, then falls back to edge_mode=’constant’.

  • ‘constant’: Return a constant specified by the constant keyword.

  • ‘wrap’: infinitely wrap the initial range in both directions.

In order for LookupTable2D to determine the wrapping period when edge_mode=’wrap’, either the x and y grid points need to be equally spaced (in which case the x-period is inferred as len(x)*(x[1]-x[0]) and similarly for y), or the first/last row/column of f must be identical, in which case the x-period is inferred as x[-1] - x[0]. (If both conditions are satisfied (equally-spaced x and y and identical first/last row/column of f, then the x-period is inferred as len(x)*(x[1]-x[0])):

>>> x = np.arange(5)
>>> y = np.arange(8)
>>> z = x + y[:, np.newaxis]  # function is x + y, dimensions of z are (8, 5)
>>> tab2d = galsim.LookupTable2D(x, y, z, edge_mode='raise')
>>> tab2d(7, 7)
ValueError: Extrapolating beyond input range.

>>> tab2d = galsim.LookupTable2D(x, y, z, edge_mode='constant', constant=1.0)
1.0

>>> tab2d = galsim.LookupTable2D(x, y, z, edge_mode='wrap')
ValueError: Cannot wrap `f` array with unequal first/last column/row.

We extend the x and y arrays with a uniform spacing, though any monotonic spacing would work. Note that the [(0,1), (0,1)] argument in np.pad below extends the z array by 0 rows/columns in the leading direction, and 1 row/column in the trailing direction:

>>> x = np.append(x, x[-1] + (x[-1]-x[-2]))
>>> y = np.append(y, y[-1] + (y[-1]-y[-2]))
>>> z = np.pad(z, [(0,1), (0,1)], mode='wrap')
>>> tab2d = galsim.LookupTable2D(x, y, z, edge_mode='wrap')
>>> tab2d(2., 2.)
4.0
>>> tab2d(2.+5, 2.)  # The period is 5 in the x direction
4.0
>>> tab2d(2.+3*5, 2.+4*8)  # The period is 8 in the y direction
4.0
Parameters:

  • x – Strictly increasing array of x positions at which to create table.

  • y – Strictly increasing array of y positions at which to create table.

  • f – Nx by Ny input array of function values.

  • dfdx – Optional first derivative of f wrt x. Only used if interpolant=’spline’. [default: None]

  • dfdy – Optional first derivative of f wrt y. Only used if interpolant=’spline’. [default: None]

  • d2fdxdy – Optional cross derivative of f wrt x and y. Only used if interpolant=’spline’. [default: None]

  • interpolant – Type of interpolation to use. One of ‘floor’, ‘ceil’, ‘nearest’, ‘linear’, ‘spline’, or a galsim.Interpolant or string convertible to one. [default: ‘linear’]

  • edge_mode – Keyword controlling how extrapolation beyond the input range is handled. See above for details. [default: ‘raise’]

  • constant – A constant to return when extrapolating beyond the input range and edge_mode='constant'. [default: 0]

__call__(x, y, grid=False)[source]

Interpolate at an arbitrary point or points.

Parameters:

  • x – Either a single x value or an array of x values at which to interpolate.

  • y – Either a single y value or an array of y values at which to interpolate.

  • grid – Optional boolean indicating that output should be a 2D array corresponding to the outer product of input values. If False (default), then the output array will be congruent to x and y.

Returns:

a scalar value if x and y are scalar, or a numpy array if x and y are arrays.

gradient(x, y, grid=False)[source]

Calculate the gradient of the function at an arbitrary point or points.

Parameters:

  • x – Either a single x value or an array of x values at which to compute the gradient.

  • y – Either a single y value or an array of y values at which to compute the gradient.

  • grid – Optional boolean indicating that output should be a 2-tuple of 2D arrays corresponding to the outer product of input values. If False (default), then the output arrays will be congruent to x and y.

Returns:

A tuple of (dfdx, dfdy) where dfdx, dfdy are single values (if x,y were single values) or numpy arrays.

galsim.table._LookupTable2D(x, y, f, interpolant, edge_mode, constant, dfdx=None, dfdy=None, d2fdxdy=None, x0=None, y0=None, xperiod=None, yperiod=None)[source]

Make a LookupTable2D but without using any of the sanity checks or array manipulation used in the normal initializer.

galsim.trapz(f, x, interpolant='linear')[source]

Integrate f(x) using the trapezoidal rule.

Equivalent to np.trapz(f,x) for 1d array inputs. Intended as a drop-in replacement, which is usually faster.

Parameters:

  • f – The ordinates of the function to integrate.

  • x – The abscissae of the function to integrate.

  • interpolant – The interpolant to use between the tabulated points. [default: ‘linear’, which matches the behavior of np.trapz]

Returns:

Estimate of the integral.